###install required packages### 
#delete hash symbol on the following lines and run (only needs to be done first time, not every time you run the code afterwards!)
#install.packages("vegan")
#install.packages("dplyr")
#install.packages("corrplot")
#install.packages("cocorresp")
#install.packages("tibble")
#install.packages("labdsv")
#install.packages("knitr")
#install.packages("dendextend")
#install.packages("pls")
#install.packages("adespatial")

Exercise 1. Introduction to R Notebooks and Markdown.

This exercise is only to get used to working with R Notebook and R Markdown. The data in this first are not multivariate and the statistical analysis is a simple linear regression. All following exercises will start with this kind of introduction.

For each exercise there will be a number of explicit questions written in green.

To pass the course, you should give your answers to these questions to one of the teachers during the computer labs.

Write your answers to the questions and any other comments here below.

Replace this text in for each question with your answers and comments…


A quick note on R Notebooks and R Markdown

These terms are often used interchangeably to mean a combination of R code, results, and formattable text, such as this one. They allow you to present your results, the code used to generate the results and your interpretation and comments all together in one document that you can easily update, and which is easy for a reader to understand (R Notebooks are a minor update of R Markdown files that allow for easier editing as you can run just part of your code and see the results in a preview. Both notebook and markdown are .Rmd files, and are written in exactly the same way, so don’t worry about the difference for this course!).


Example of a simple analysis

First we load or create the data we will be using. The box below is an example of a code chunk.

df <- read.table(header = TRUE, sep = ",", 
               text = "x,y 
               1,2
               1.5,3.5
               2,5
               2.5,5
               3,7")

# Anything you write in a code chunk after a hash mark becomes a comment like this
# and is ignored when running the code. If you are working in an R markdown notebook
# such as this however, you can also write your longer notes in the text sections 
# outside the code chunks.

So I could write my notes here instead, which is probably easier to read for anything more than a few words of explanation.

Now we can do a simple linear regression. The results are stored in the object “reg”.

reg <- lm(y~x, df) 

Plot the result.

plot(y~x, df)
lines(df$x, predict(reg), col = 'red') # Add a regression line to the plot

Get a summary of the results

summary(reg)
## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Residuals:
##     1     2     3     4     5 
## -0.20  0.15  0.50 -0.65  0.20 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -0.1000     0.6819  -0.147  0.89271   
## x             2.3000     0.3215   7.155  0.00562 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5083 on 3 degrees of freedom
## Multiple R-squared:  0.9446, Adjusted R-squared:  0.9262 
## F-statistic: 51.19 on 1 and 3 DF,  p-value: 0.005622

The R packages and data we will be using.

The “dune” dataset we will use is included in the vegan package (an extensive collection of tools for multivariate analysis), but we will also be using some slightly modified and extended versions of it, which we will load as needed. We will also load some other packages to extend the range of analyses available to us.




##Q2

Question 2: Simple exercise for illustrating effects of different distance metrics. Should not be based on ordination or classification.

Simple exercise for illustrating effects of different distance metrics. Should not be based on ordination or classification.

#TBC, Martyn




##Q3

Question 3: In this exercise we will use Principal Components Analysis on the dune meadow explanatory data. We will also explore the effects of different kinds of scaling on the ordination plot. NOTE that the scaling issues applies to all ordination techniques!

Compare the plots on the PCA with and without scaling (PCA plot 4) and discuss the differences. Is there any way to see that one is on scaled data and one is not? What is the difference between PCA plot 4 and PCA plot 5? Which one is best to use? Discuss the difference between the plots in plot PCA plot 6 and 7!

library(tidyverse)
library(vegan)
library(corrplot)

data("dune")
data("dune.env")
dummy_management <- as.data.frame(model.matrix( ~ Management - 1, data=dune.env )) 
#add these to the dataset
dune.env.original <- dune.env #we keep a copy of the original version
dune.env <- dune.env %>% select(A1, Moisture, Manure, Use) %>% cbind(.,dummy_management) 
dune.env$Moisture <- as.numeric(as.character(dune.env$Moisture)) #make numeric
dune.env$Manure <- as.numeric(as.character(dune.env$Manure))
dune.env$Use <- as.numeric(dune.env$Use)
#make column names shorter
dune.env <- dune.env %>% rename(BF = ManagementBF, HF = ManagementHF,
                                NM = ManagementNM, SF = ManagementSF)

## PCA plot 1:
# Start by looking at pairwise correlations among the variables.
dune.env_cor<-cor(dune.env, method = "kendall")
corrplot(dune.env_cor, type = "upper", order = "hclust", tl.col = "black", tl.srt = 45)

## PCA plot 2:
# Then make a PCA, with scaling of the data
env.pca <- rda(dune.env, scale = TRUE) #vegan uses the same function for PCA and RDA, just depends on if it is constrained or not.
biplot(env.pca) #plot the results using the default plot scaling which is "species"

env.pca #summarise results
## Call: rda(X = dune.env, scale = TRUE)
## 
##               Inertia Rank
## Total               8     
## Unconstrained       8    7
## Inertia is correlations 
## 
## Eigenvalues for unconstrained axes:
##    PC1    PC2    PC3    PC4    PC5    PC6    PC7 
## 2.7348 1.9195 1.3056 1.0477 0.5106 0.4141 0.0677
summary(eigenvals(env.pca)) #see variance explained
## Importance of components:
##                          PC1    PC2    PC3   PC4     PC5     PC6      PC7
## Eigenvalue            2.7348 1.9195 1.3056 1.048 0.51059 0.41413 0.067668
## Proportion Explained  0.3418 0.2399 0.1632 0.131 0.06382 0.05177 0.008458
## Cumulative Proportion 0.3418 0.5818 0.7450 0.876 0.93977 0.99154 1.000000
## PCA plot 3:
# Continue with a PCA on the same data, but without scaling of the data
env.pca2 <- rda(dune.env, scale = FALSE) #vegan uses the same function for PCA and RDA, just depends on if it is constrained or not.
biplot(env.pca2)#plot the results

env.pca2 #summarise results
## Call: rda(X = dune.env, scale = FALSE)
## 
##               Inertia Rank
## Total           11.59     
## Unconstrained   11.59    7
## Inertia is variance 
## 
## Eigenvalues for unconstrained axes:
##   PC1   PC2   PC3   PC4   PC5   PC6   PC7 
## 6.436 2.314 1.993 0.435 0.255 0.126 0.035
summary(eigenvals(env.pca2)) #see variance explained
## Importance of components:
##                          PC1    PC2    PC3     PC4     PC5     PC6      PC7
## Eigenvalue            6.4360 2.3137 1.9929 0.43530 0.25488 0.12593 0.034957
## Proportion Explained  0.5551 0.1996 0.1719 0.03755 0.02198 0.01086 0.003015
## Cumulative Proportion 0.5551 0.7547 0.9266 0.96414 0.98612 0.99698 1.000000
## PCA plot 4:
# Compare the PCAs with and without scaling of data
# Default scaling of the plots, which is "species"
par(mfrow = c(1,2))
biplot(env.pca, main = "PCA with scaling of data")
biplot(env.pca2, main = "PCA without scaling of data")

par(mfrow = c(1,1))

## PCA plot 5:
# Plots from the same PCA, but with plot scaling focused on sites and on species
par(mfrow = c(1,2))
biplot(env.pca, scaling = "sites", main = "Plot scaling on sites")
biplot(env.pca, scaling = "species", main = "Plot scaling on species")

par(mfrow = c(1,1))

## PCA plot 6:
# Then, a plot from the PCA on un-scaled data, 
# but with plot scaling to mimic a PCA on scaled data (correlation = TRUE)
# Second plot is the PCA on scaled data
# Both plots with plot scaling focused on species
par(mfrow = c(1,2))
biplot(env.pca2, scaling = "species", correlation = TRUE, main = "PCA on un-sclaed data,\nplot scaling on correlations")
biplot(env.pca, scaling = "species", correlation = FALSE, main = "PCA on scaled data")

par(mfrow = c(1,1))

## PCA plot 7:
# Finally a plot from the PCA on un-scaled data,
# but with plot scaling to mimic a PCA on scaled data
# Second plot is the PCA on scaled data
# Both plots with plot scaling focused on sites
par(mfrow = c(1,2))
biplot(env.pca2, scaling = "sites", correlation = TRUE, main = "PCA on un-scaled data,\nplot scaling on correlations")
biplot(env.pca, scaling = "sites", main = "PCA on scaled data")

par(mfrow = c(1,1))




##Q4

Question 4: Do a CA-ordination on the Dune Meadow species dataset.What are the results telling you?

Give a conceptual description on why objects/samples and descriptors/species to the left differ from objects and descriptors to the right, and those at the bottom from those at the top! (Plant ecologists may give a more detailed description, using their knowledge about the species in the dataset).

dune.ca <- cca(dune)
plot(dune.ca)

dune.ca
## Call: cca(X = dune)
## 
##               Inertia Rank
## Total           2.115     
## Unconstrained   2.115   19
## Inertia is scaled Chi-square 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.5360 0.4001 0.2598 0.1760 0.1448 0.1079 0.0925 0.0809 
## (Showing 8 of 19 unconstrained eigenvalues)
summary(eigenvals(dune.ca)) #proportion variance explained
## Importance of components:
##                          CA1    CA2    CA3     CA4     CA5     CA6     CA7
## Eigenvalue            0.5360 0.4001 0.2598 0.17598 0.14476 0.10791 0.09247
## Proportion Explained  0.2534 0.1892 0.1228 0.08319 0.06844 0.05102 0.04372
## Cumulative Proportion 0.2534 0.4426 0.5654 0.64858 0.71702 0.76804 0.81175
##                           CA8     CA9    CA10    CA11    CA12    CA13     CA14
## Eigenvalue            0.08091 0.07332 0.05630 0.04826 0.04125 0.03523 0.020529
## Proportion Explained  0.03825 0.03466 0.02661 0.02282 0.01950 0.01665 0.009705
## Cumulative Proportion 0.85000 0.88467 0.91128 0.93410 0.95360 0.97025 0.979955
##                           CA15     CA16     CA17     CA18     CA19
## Eigenvalue            0.014911 0.009074 0.007938 0.007002 0.003477
## Proportion Explained  0.007049 0.004290 0.003753 0.003310 0.001644
## Cumulative Proportion 0.987004 0.991293 0.995046 0.998356 1.000000




##Q5

Question 5: Repeat exercise 4, but with DCA ordination instead.

Look at the eigenvalues, the length of gradient, the total variation and the ordination diagram. Explain the differences between results from CA and DCA.

dune.dca <- decorana(dune)
dune.dca 
## 
## Call:
## decorana(veg = dune) 
## 
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
## 
##                   DCA1   DCA2    DCA3    DCA4
## Eigenvalues     0.5117 0.3036 0.12125 0.14267
## Decorana values 0.5360 0.2869 0.08136 0.04814
## Axis lengths    3.7004 3.1166 1.30055 1.47888
#Detrended correspondence analysis (function decorana).
#Note that we do not get "variation explained" in the R implementation of DCA (and some other functions).
#Here, the developer explains why, "The total amount of variation is undefined in detrended 
#correspondence analysis and therefore proportions from total are unknown and undefined. 
#DCA is not a method for decomposition of variation, and therefore
#these proportions would not make sense either.
plot(dune.dca)




##Q6

Question 6: Add one new species with cover degree 9 at the site with the lowest number of species. Do a new DCA-ordination.

What happens with the eigenvalues and ordination diagram?

# calculate the number of species present in each plot
rowSums(dune != 0)
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 
##  5 10 10 13 14 11 13 12 13 12  9  9 10  7  8  8  7  9  9  8
# find the plot with the lowest number of species
which.min(rowSums(dune != 0))
## 1 
## 1
# create a new data frame with a new species with cover degree 9 at site 1
NewSp <- c(9, rep(0,19))
SpeciesNewSp <- data.frame(dune, NewSp)
# run a new DCA with the extra species
DCANewSp <- decorana(SpeciesNewSp)
DCANewSp
## 
## Call:
## decorana(veg = SpeciesNewSp) 
## 
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
## 
##                   DCA1   DCA2    DCA3    DCA4
## Eigenvalues     0.5229 0.3352 0.13870 0.09294
## Decorana values 0.5436 0.2636 0.08678 0.04037
## Axis lengths    3.4890 3.6245 1.49207 1.24554
# plot the DCA
ordiplot(DCANewSp, type="text", main="DCA with new sp.")

You may be surprised that you haven’t got any total inertia values when printing decorana results, although in other software (e.g. CANOCO) these are available, together with percentage variance explained by particular axes. The reason for this is that DCA does not support the concept of total inertia values (also, it produces only four axes, i.e. four eigenvalues). See the comment in the code for the question 5.




##Q7

### Question 7: Effects from increased variation. Add one new sample to the data, then add a cover of 9 for the existing species Chenopodium album “Che_alb” in the new sample. Make a new DCA-ordination.

Next, remove the Che alb = 9, and add instead a cover of 5 for the exisiting species Leo_aut and make a new DCA.

What happens with the eigenvalues and ordination after adding Che alb and after adding Leo aut? Why did we choose to add these two species in the new sample?

# figure out the column number for Chenopodium album (abbreviated "Chenalbu")
which(colnames(dune)=="Che_alb")
## integer(0)
# create a new data frame with a new sample containing Chenopodium album with a cover of 9
NewSample <- c(rep(0,7), 9, rep(0,22))
SpeciesNewSample <- data.frame(rbind(dune, NewSample))
# run a new DCA with the extra sample
DCANewSample <- decorana(SpeciesNewSample)
DCANewSample
## 
## Call:
## decorana(veg = SpeciesNewSample) 
## 
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
## 
##                   DCA1   DCA2   DCA3    DCA4
## Eigenvalues     0.9012 0.5109 0.3028 0.13538
## Decorana values 0.9020 0.5345 0.2856 0.09427
## Axis lengths    1.7784 3.6794 3.1985 1.51095
# plot the DCA
ordiplot(DCANewSample, type="text", main="DCA with new sample Chenopodium album")

# figure out the column number for Leontodon autumnalis (abbreviated "Leo aut")
which(colnames(dune)=="Leo_aut")
## integer(0)
# create a new data frame with a new sample containing Leontodon autumnalis with a cover of 5
NewSample <- c(rep(0,15), 5, rep(0,10))
SpeciesNewSample <- data.frame(rbind(dune, NewSample))
# run a new DCA with the extra sample
DCANewSample <- decorana(SpeciesNewSample)
DCANewSample
## 
## Call:
## decorana(veg = SpeciesNewSample) 
## 
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
## 
##                   DCA1   DCA2   DCA3    DCA4
## Eigenvalues     0.5125 0.3040 0.3030 0.12779
## Decorana values 0.5361 0.2869 0.2187 0.07207
## Axis lengths    3.7039 3.0498 0.5132 1.44059
# plot the DCA
ordiplot(DCANewSample, type="text", main="DCA with new sample Leontodon autumnalis")




##Q8

Question 8: Non Metric Multidimensional Scaling

Compare the NMDS ordination with the ordinations obtained using CA and DCA! What is the stress value for the NMDS?

# NMDS on original data set with 2 reduced dimensions
NMDS <- metaMDS(dune, k=2)
## Run 0 stress 0.1192678 
## Run 1 stress 0.1812936 
## Run 2 stress 0.1192678 
## ... New best solution
## ... Procrustes: rmse 3.554456e-05  max resid 0.0001037297 
## ... Similar to previous best
## Run 3 stress 0.2341478 
## Run 4 stress 0.1808914 
## Run 5 stress 0.180892 
## Run 6 stress 0.1192678 
## ... New best solution
## ... Procrustes: rmse 2.856648e-05  max resid 8.100342e-05 
## ... Similar to previous best
## Run 7 stress 0.1192678 
## ... Procrustes: rmse 4.60395e-05  max resid 0.0001366484 
## ... Similar to previous best
## Run 8 stress 0.1808911 
## Run 9 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 0.02027803  max resid 0.06500843 
## Run 10 stress 0.1886532 
## Run 11 stress 0.119269 
## Run 12 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 2.017647e-05  max resid 6.602145e-05 
## ... Similar to previous best
## Run 13 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 7.419144e-06  max resid 2.567168e-05 
## ... Similar to previous best
## Run 14 stress 0.119268 
## Run 15 stress 0.1183186 
## ... Procrustes: rmse 2.117589e-05  max resid 5.801426e-05 
## ... Similar to previous best
## Run 16 stress 0.1192678 
## Run 17 stress 0.1183186 
## ... Procrustes: rmse 2.214182e-05  max resid 6.882882e-05 
## ... Similar to previous best
## Run 18 stress 0.119268 
## Run 19 stress 0.119268 
## Run 20 stress 0.1183186 
## ... New best solution
## ... Procrustes: rmse 2.785698e-06  max resid 8.22521e-06 
## ... Similar to previous best
## *** Solution reached
NMDS
## 
## Call:
## metaMDS(comm = dune, k = 2) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     dune 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1183186 
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'dune'
# plot the NMDS
ordiplot(NMDS, type="text", main="NMDS")

#plot the CA
ordiplot(dune.ca, type="text", main="DCA")

#plot the DCA
ordiplot(dune.dca, type="text", main="DCA")



##Q9

Question 9: PCoA

data quality


#To follow, Martyn

##Q10

Question 10: Data quality assessment

data quality

# To follow, Ulf




##Q11

Question 11: Do an RDA and a CCA on the Dune Meadow data set, using the full set of environmental variables.

Compare the ordination diagrams and numerical results! What does the significance test tell you? When only looking at the results, which method explains the variation best? Taking the nature of the data into account, which method should be used?

# RDA 
dune.rda <- rda(dune ~ . , dune.env.original)
summary(dune.rda)
## 
## Call:
## rda(formula = dune ~ A1 + Moisture + Management + Use + Manure,      data = dune.env.original) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total           84.12     1.0000
## Constrained     63.21     0.7513
## Unconstrained   20.92     0.2487
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                          RDA1    RDA2    RDA3   RDA4   RDA5    RDA6    RDA7
## Eigenvalue            22.3955 16.2076 7.03891 4.0380 3.7602 2.60874 2.16693
## Proportion Explained   0.2662  0.1927 0.08367 0.0480 0.0447 0.03101 0.02576
## Cumulative Proportion  0.2662  0.4589 0.54256 0.5906 0.6353 0.66627 0.69203
##                          RDA8    RDA9   RDA10    RDA11    RDA12     PC1     PC2
## Eigenvalue            1.80327 1.40421 0.91739 0.581545 0.283927 6.62687 4.30914
## Proportion Explained  0.02144 0.01669 0.01091 0.006913 0.003375 0.07878 0.05122
## Cumulative Proportion 0.71346 0.73015 0.74106 0.747973 0.751348 0.83012 0.88135
##                           PC3     PC4     PC5    PC6      PC7
## Eigenvalue            3.54913 2.54648 2.34027 0.9335 0.612096
## Proportion Explained  0.04219 0.03027 0.02782 0.0111 0.007276
## Cumulative Proportion 0.92354 0.95381 0.98163 0.9927 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                          RDA1    RDA2   RDA3    RDA4    RDA5    RDA6    RDA7
## Eigenvalue            22.3955 16.2076 7.0389 4.03795 3.76017 2.60874 2.16693
## Proportion Explained   0.3543  0.2564 0.1114 0.06389 0.05949 0.04127 0.03428
## Cumulative Proportion  0.3543  0.6107 0.7221 0.78600 0.84549 0.88676 0.92105
##                          RDA8    RDA9   RDA10    RDA11    RDA12
## Eigenvalue            1.80327 1.40421 0.91739 0.581545 0.283927
## Proportion Explained  0.02853 0.02222 0.01451 0.009201 0.004492
## Cumulative Proportion 0.94958 0.97179 0.98631 0.995508 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  6.322924 
## 
## 
## Species scores
## 
##              RDA1     RDA2      RDA3      RDA4      RDA5      RDA6
## Achimill  0.58181  0.08704 -0.105553 -0.043694 -0.422290 -0.062221
## Agrostol -1.16975 -0.94792  0.060309  0.334142  0.030721  0.222829
## Airaprae -0.11716  0.20786  0.018154 -0.180185 -0.050376  0.041993
## Alopgeni -0.50957 -1.41677 -0.322262 -0.116416 -0.013318 -0.300469
## Anthodor  0.39549  0.45832 -0.411832 -0.136384 -0.336365  0.263272
## Bellpere  0.36549 -0.10133  0.162631 -0.160330 -0.370665  0.137716
## Bromhord  0.53405 -0.19418  0.161217 -0.048422 -0.585999 -0.030597
## Chenalbu -0.03146 -0.03881 -0.034281  0.004236 -0.017102  0.003806
## Cirsarve  0.02198 -0.09707  0.075983  0.009750  0.017309  0.069692
## Comapalu -0.19434  0.09464  0.045000  0.157251 -0.021772 -0.021112
## Eleopalu -1.07946  0.08439  0.006165  0.698662 -0.258796  0.127090
## Elymrepe  0.57194 -0.67735  0.299352 -0.407119  0.233601  0.408690
## Empenigr -0.07889  0.08081  0.032442 -0.076400 -0.005971  0.030773
## Hyporadi -0.12994  0.36832  0.135169 -0.148462  0.009242 -0.069277
## Juncarti -0.46506 -0.11361 -0.043146  0.023106  0.051743  0.032617
## Juncbufo -0.13148 -0.45491 -0.399997 -0.320717  0.288602 -0.369510
## Lolipere  1.57085 -0.39016  0.727211  0.554451  0.340886 -0.082511
## Planlanc  0.89986  0.57656 -0.564648  0.273286  0.038960  0.120162
## Poaprat   0.98798 -0.44183  0.418705  0.121187  0.153703  0.036591
## Poatriv   0.74253 -1.39782 -0.559135  0.008607 -0.428312  0.162203
## Ranuflam -0.55236  0.08108  0.004173  0.185963 -0.131421  0.083262
## Rumeacet  0.57938  0.06076 -0.940669  0.137799  0.309740  0.163586
## Sagiproc -0.21299 -0.45646  0.011857 -0.174011  0.208064 -0.216092
## Salirepe -0.28603  0.50726  0.148922 -0.456522  0.030091  0.164904
## Scorautu  0.33487  0.50344  0.086424 -0.148075 -0.258282 -0.288324
## Trifprat  0.42084  0.16346 -0.528815  0.260068  0.196986  0.120890
## Trifrepe  0.39232  0.08889 -0.146355  0.260759 -0.280784 -0.613465
## Vicilath  0.11379  0.13700  0.111709  0.049534  0.011107 -0.146626
## Bracruta -0.01603  0.41968 -0.395909  0.008729  0.367890 -0.142926
## Callcusp -0.52015  0.21794  0.061946  0.231494 -0.019057  0.025981
## 
## 
## Site scores (weighted sums of species scores)
## 
##        RDA1    RDA2     RDA3     RDA4     RDA5    RDA6
## 1   1.01596 -0.1799  2.59463 -0.11443  2.05964  2.1377
## 2   1.72853 -1.0870  1.31739 -0.89469 -3.45007 -1.0625
## 3   0.61315 -2.4671  1.08784 -0.06821  0.64265  0.3666
## 4   0.12048 -2.1589  1.71771 -0.07977  0.48648  1.8652
## 5   1.68284  0.4191 -2.55443 -0.80052 -1.22481  3.3283
## 6   2.08298  1.2151 -3.01397  2.45080  1.69576 -0.3949
## 7   1.89757  0.4084 -1.15310  1.15696  0.19469  0.1648
## 8  -0.79508 -1.2469  0.59088  1.85859  0.03913 -0.7791
## 9   0.08193 -1.7029 -0.41962 -2.00241  1.90642  0.0328
## 10  1.95184  0.6210  0.49322  0.87955 -3.21475 -1.7550
## 11  1.07928  1.3401  1.77536  0.93077  2.20291 -2.9068
## 12 -1.19429 -1.6563 -2.36651 -1.91407  1.14492 -3.7773
## 13 -0.82017 -2.3013 -1.46058 -1.16242 -1.40595 -0.6812
## 14 -1.69462  1.0179  0.55226  2.12686 -1.38338 -1.1655
## 15 -1.92197  1.0316  0.18669  1.38332  0.21043  0.9451
## 16 -2.81781 -0.4853 -0.45260  2.98066 -0.24307  2.1409
## 17 -0.26124  1.7408 -0.04449 -2.00579 -0.85795  1.4187
## 18  0.25044  1.9245  0.47979 -1.27429  1.11822 -1.1395
## 19 -0.64055  2.2851  0.18205 -3.72278 -0.40927 -1.1779
## 20 -2.35927  1.2820  0.48747  0.27188  0.48801  2.4397
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##        RDA1     RDA2     RDA3     RDA4     RDA5     RDA6
## 1   0.87487 -0.90608  2.36103 -0.66159  2.53531  1.55448
## 2   1.91194 -0.62145  1.07037 -0.88335 -3.00393 -2.04867
## 3   0.43540 -1.95933  1.52026  0.20428  0.30748  1.42187
## 4   0.43933 -1.94048  1.51888  0.19490  0.34600  1.39312
## 5   1.65691  0.51332 -2.47861  0.06329 -1.37887  3.19410
## 6   2.10890  1.12091 -3.08978  1.58700  1.84982 -0.26074
## 7   1.48330 -0.04811 -0.36778  1.16790  0.69199 -0.12568
## 8  -0.56422 -1.25593  0.05259  1.83632 -0.90431  0.49765
## 9   0.26534 -1.23731 -0.66664 -1.99108  2.35256 -0.95341
## 10  1.52623 -0.09659  0.82501 -0.36238 -2.76909 -0.83620
## 11  1.32148  1.59213  1.69059  2.16136  1.31111 -2.83941
## 12 -1.37770 -2.12186 -2.11948 -1.92541  0.69878 -2.79106
## 13 -1.25782 -1.55141 -1.37052  0.16936 -0.68373  0.15215
## 14 -2.37447  1.51577  0.72709  2.22806  0.23064 -0.39953
## 15 -1.51025  0.37596  0.17244  0.91534 -0.66586 -0.02248
## 16 -2.19675 -0.76968 -0.78968  1.66022 -0.51916  0.32129
## 17  0.02328  1.73218 -0.60987 -1.31102 -0.82795 -0.08329
## 18  0.38018  2.38966  0.25985 -1.97999  0.59091  0.65301
## 19 -1.57691  1.61530  0.64851 -1.52722 -0.11937  0.61515
## 20 -1.56904  1.65300  0.64575 -1.54597 -0.04232  0.55766
## 
## 
## Biplot scores for constraining variables
## 
##                 RDA1     RDA2      RDA3     RDA4      RDA5     RDA6
## A1           -0.5361  0.06882 -0.218452  0.31811 -0.146216  0.05232
## Moisture.L   -0.9078 -0.11664 -0.056742  0.05789 -0.136471 -0.01179
## Moisture.Q   -0.1012  0.43397 -0.036128  0.40058 -0.008317  0.14269
## Moisture.C   -0.1978  0.01609  0.364691  0.19843 -0.452618  0.35989
## ManagementHF  0.4043 -0.07409 -0.534962  0.21752  0.213258  0.19208
## ManagementNM -0.5114  0.71630  0.142286 -0.24855 -0.064357  0.10191
## ManagementSF -0.2379 -0.71375  0.086470 -0.02765  0.207180  0.15834
## Use.L        -0.1398 -0.25485  0.123795  0.81530  0.138325 -0.26210
## Use.Q        -0.1173  0.54889 -0.004563  0.02884 -0.099898  0.04341
## Manure.L      0.2718 -0.78470  0.068253  0.26784  0.140391  0.28325
## Manure.Q     -0.3913  0.30047  0.585951 -0.19204  0.183926  0.39222
## Manure.C      0.5204 -0.16720  0.322839 -0.18763  0.228598 -0.20896
## 
## 
## Centroids for factor constraints
## 
##                 RDA1     RDA2      RDA3     RDA4    RDA5     RDA6
## Moisture1     1.3911  0.57720 -0.079190  0.20780  0.3709  0.01816
## Moisture2     0.6061 -0.56605  0.813569 -0.31856 -0.7359  0.47388
## Moisture4    -0.5562 -1.67959 -1.393063 -1.95824  1.5257 -1.87223
## Moisture5    -1.5785  0.22614  0.012311  0.53373 -0.3863  0.24598
## ManagementBF  1.5865  0.29136  1.195323  0.30521 -1.4873 -1.90809
## ManagementHF  0.9900 -0.18143 -1.310046  0.53268  0.5222  0.47038
## ManagementNM -1.1045  1.54698  0.307294 -0.53680 -0.1390  0.22008
## ManagementSF -0.5138 -1.54147  0.186749 -0.05971  0.4474  0.34198
## UseHayfield   0.1009  0.93851 -0.196574 -1.23634 -0.3134  0.44957
## UseHaypastu   0.2031 -0.95047  0.007901 -0.04993  0.1730 -0.07517
## UsePasture   -0.4661  0.20684  0.262562  1.81077  0.1621 -0.50914
## Manure1       1.0377  0.08608  0.616318 -0.06403  0.2982 -1.54300
## Manure2       1.0750 -0.27727 -1.654376 -0.28962 -0.4585 -0.47659
## Manure3      -0.6339 -0.90628 -0.618849  1.20845 -0.3538  0.21135
## Manure4       0.5832 -1.60196  1.800060 -0.08747  1.0629  1.45649
anova(dune.rda)
# CCA 
dune.cca <- cca(dune ~ . , dune.env.original)
summary(dune.cca)
## 
## Call:
## cca(formula = dune ~ A1 + Moisture + Management + Use + Manure,      data = dune.env.original) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.1153     1.0000
## Constrained    1.5032     0.7106
## Unconstrained  0.6121     0.2894
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                         CCA1   CCA2    CCA3    CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4671 0.3410 0.17606 0.15317 0.09528 0.07027 0.05887
## Proportion Explained  0.2208 0.1612 0.08323 0.07241 0.04504 0.03322 0.02783
## Cumulative Proportion 0.2208 0.3821 0.46529 0.53770 0.58275 0.61596 0.64379
##                          CCA8    CCA9   CCA10   CCA11    CCA12    CA1     CA2
## Eigenvalue            0.04993 0.03183 0.02596 0.02282 0.010818 0.2724 0.10876
## Proportion Explained  0.02360 0.01505 0.01227 0.01079 0.005114 0.1288 0.05142
## Cumulative Proportion 0.66740 0.68245 0.69472 0.70551 0.710627 0.8394 0.89081
##                           CA3     CA4     CA5     CA6      CA7
## Eigenvalue            0.08975 0.06305 0.03489 0.02529 0.017984
## Proportion Explained  0.04243 0.02981 0.01649 0.01196 0.008502
## Cumulative Proportion 0.93324 0.96305 0.97954 0.99150 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1   CCA2   CCA3   CCA4    CCA5    CCA6    CCA7
## Eigenvalue            0.4671 0.3410 0.1761 0.1532 0.09528 0.07027 0.05887
## Proportion Explained  0.3108 0.2269 0.1171 0.1019 0.06339 0.04674 0.03916
## Cumulative Proportion 0.3108 0.5376 0.6548 0.7567 0.82005 0.86679 0.90595
##                          CCA8    CCA9   CCA10   CCA11    CCA12
## Eigenvalue            0.04993 0.03183 0.02596 0.02282 0.010818
## Proportion Explained  0.03322 0.02118 0.01727 0.01518 0.007197
## Cumulative Proportion 0.93917 0.96035 0.97762 0.99280 1.000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##               CCA1     CCA2      CCA3      CCA4      CCA5      CCA6
## Achimill  0.857707  0.49484 -0.045074 -0.057950  0.670679  0.056946
## Agrostol -0.750987 -0.51041 -0.001514 -0.006522  0.110328 -0.113279
## Airaprae -0.786606  1.84744  0.777309  1.068417  0.558936 -0.142093
## Alopgeni -0.296764 -0.98554  0.015875  0.454024  0.295335  0.124614
## Anthodor  0.366346  1.04925 -0.217724  0.412435  0.563903 -0.226744
## Bellpere  0.707119  0.06620  0.500239 -0.369709  0.592272 -0.330213
## Bromhord  0.853359 -0.02616  0.335626 -0.427024  0.865565 -0.093885
## Chenalbu -0.672293 -1.27568 -0.417041  0.994310  0.640906  0.521336
## Cirsarve  0.386092 -0.93251  0.963978 -0.468823  0.239146 -0.908823
## Comapalu -2.015388  0.10446 -1.185954 -1.972940 -0.083546 -0.219456
## Eleopalu -1.515611 -0.03069 -0.455387 -0.508862  0.111183 -0.132850
## Elymrepe  0.610453 -0.66152  0.628614 -0.124262 -0.491387 -0.510855
## Empenigr -1.378845  1.51426  1.387645  0.973499 -0.434196 -0.642379
## Hyporadi -0.655401  1.49855  0.738512  0.362909 -0.104564  0.310460
## Juncarti -0.948447 -0.19092 -0.028958  0.148431 -0.292973 -0.155408
## Juncbufo -0.064392 -1.01588 -0.297937  1.118684 -0.196979  0.823725
## Lolipere  0.607215 -0.09885  0.262228 -0.409402 -0.312573  0.082660
## Planlanc  0.800409  0.77629 -0.606615 -0.038048 -0.021902  0.152465
## Poaprat   0.471703 -0.13797  0.252857 -0.205726 -0.176841  0.061944
## Poatriv   0.303929 -0.50098 -0.033488  0.157680  0.180752 -0.110423
## Ranuflam -1.384091  0.11712 -0.129524 -0.120128  0.041536 -0.207804
## Rumeacet  0.791036  0.13052 -1.100806  0.537567 -0.416656 -0.396323
## Sagiproc -0.290759 -0.52139  0.288191  0.373446 -0.109281  0.193977
## Salirepe -0.889942  1.59655  1.182347  0.817665 -0.457263 -0.271421
## Scorautu -0.009838  0.37963  0.098333 -0.075560  0.001781  0.125526
## Trifprat  1.054359  0.56735 -1.455764  0.398295 -0.539708 -0.564846
## Trifrepe  0.014698  0.11386 -0.184205 -0.244273  0.089510  0.207091
## Vicilath  0.511432  0.88591  0.273760 -1.050861 -0.335122  1.591159
## Bracruta -0.147674  0.27946 -0.137684  0.084272 -0.210685  0.147748
## Callcusp -1.773372  0.55553 -0.188876 -0.426816 -0.097141 -0.008786
## 
## 
## Site scores (weighted averages of species scores)
## 
##        CCA1     CCA2      CCA3    CCA4     CCA5     CCA6
## 1   1.19460 -0.71633  1.656429 -1.4249 -2.23242 -1.09183
## 2   0.86805 -0.35777  0.934868 -0.9242  1.33133 -0.48738
## 3   0.33794 -1.04078  0.824936 -0.1786 -0.04887 -0.42012
## 4   0.25095 -0.99294  1.237259 -0.4411  0.35593 -1.14729
## 5   1.11991  0.45932 -1.022763  0.2895  0.27912 -1.76102
## 6   0.99305  0.73388 -2.001441  0.3758 -1.04277 -0.73162
## 7   1.03098  0.34363 -1.083231  0.2079 -0.14394  0.29945
## 8  -0.66671 -0.71037  0.004385 -0.1742 -0.10923  0.07552
## 9   0.09269 -1.09341  0.218954  0.9201 -1.32966  0.09022
## 10  0.95315  0.58996  0.146744 -0.9667  1.81728  0.84400
## 11  0.47318  0.74856  0.535395 -1.0620 -1.51667  3.25324
## 12 -0.27934 -1.30695 -0.512852  2.0022  0.26928  1.96951
## 13 -0.37400 -1.45815 -0.153034  1.4764  1.08084  1.21481
## 14 -2.04173  0.23744 -1.448971 -2.6038  0.41722 -0.22568
## 15 -1.93799 -0.04255 -1.367867 -1.8081 -0.32644 -0.85946
## 16 -1.91272 -0.56130 -0.881308 -0.5388  0.40537 -0.81739
## 17  0.33533  2.74717  0.494859  1.7252  3.00069  0.15330
## 18  0.26805  1.23029  0.753438 -0.3018 -1.09614  1.50399
## 19 -0.75573  2.53794  2.225936  2.4221  0.13516  0.14348
## 20 -2.03905  0.80938  0.463176  0.2188 -1.15999 -1.41953
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##       CCA1    CCA2      CCA3     CCA4    CCA5     CCA6
## 1   0.7577 -1.0533  2.075818 -0.62357 -2.6193 -0.57163
## 2   1.0307 -0.2511  0.942759 -0.72795  1.0939  0.10633
## 3   0.3825 -0.9571  0.943652 -0.50806  0.2566 -0.92249
## 4   0.3861 -0.9325  0.963978 -0.46882  0.2391 -0.90882
## 5   1.0281  0.3743 -1.180481 -0.15376  0.1863 -2.04984
## 6   1.0753  0.8100 -1.860152  0.77287 -0.9596 -0.47288
## 7   1.0282  0.1537 -0.720078  0.01393 -0.2159  0.69023
## 8  -0.8347 -0.6324 -0.367053 -0.18633  0.2121 -0.93866
## 9   0.2553 -0.9868  0.226845  1.11630 -1.5671  0.68393
## 10  0.8517  0.5884  0.543136 -0.45186  1.9874  0.90302
## 11  0.3961  0.6106 -0.007613 -2.01145 -1.4337  2.39468
## 12 -0.4745 -1.4349 -0.522321  1.76672  0.5542  1.25706
## 13 -0.6723 -1.2757 -0.417041  0.99431  0.6409  0.52134
## 14 -2.3490  0.6722 -1.087969 -2.00082 -0.2983 -0.03446
## 15 -1.6817 -0.4633 -1.283939 -1.94506  0.1312 -0.40445
## 16 -1.4075 -0.6080 -0.607258  0.19306  0.5431  0.63171
## 17  0.1018  2.3472 -0.138197  1.21079  2.0486  0.60834
## 18  0.4019  1.7339  0.567131  0.27132 -0.4606  0.67225
## 19 -1.3788  1.5143  1.387645  0.97350 -0.4342 -0.64238
## 20 -1.3717  1.5635  1.428297  1.05197 -0.4691 -0.61505
## 
## 
## Biplot scores for constraining variables
## 
##                 CCA1     CCA2     CCA3     CCA4     CCA5     CCA6
## A1           -0.5534 -0.14431 -0.51439 -0.29750  0.10556 -0.18207
## Moisture.L   -0.9117 -0.18833  0.03981  0.20831  0.10483 -0.03496
## Moisture.Q   -0.2030  0.34831 -0.29281 -0.20645 -0.28387 -0.13407
## Moisture.C   -0.2100  0.10223  0.30648 -0.29962  0.52384 -0.28198
## ManagementHF  0.3604 -0.01771 -0.54879  0.22226 -0.32807 -0.29032
## ManagementNM -0.5914  0.64415  0.15809 -0.01577 -0.04877 -0.07689
## ManagementSF -0.1243 -0.67548  0.21199  0.15443  0.10327 -0.01928
## Use.L        -0.1943 -0.34773 -0.33160 -0.42673 -0.10227  0.22872
## Use.Q        -0.1889  0.51516 -0.02141 -0.01123 -0.05134  0.16126
## Manure.L      0.3295 -0.69510  0.03992 -0.02906  0.04350 -0.26625
## Manure.Q     -0.4229  0.22292  0.52833 -0.20228 -0.12773 -0.23201
## Manure.C      0.4853 -0.11991  0.31222 -0.19850 -0.15170  0.17650
## 
## 
## Centroids for factor constraints
## 
##                  CCA1     CCA2     CCA3      CCA4     CCA5    CCA6
## Moisture1     0.86961  0.39190 -0.30729 -0.271184 -0.42482  0.0229
## Moisture2     0.49529 -0.13803  0.71613 -0.298514  0.95955 -0.2087
## Moisture4    -0.07641 -1.19047 -0.11369  1.411949 -0.60287  0.9444
## Moisture5    -1.31323  0.06247 -0.07829  0.008202  0.07168 -0.2261
## ManagementBF  0.79133  0.29313  0.53596 -0.977522  0.73099  1.0250
## ManagementHF  0.53656 -0.02637 -0.81693  0.330865 -0.48837 -0.4322
## ManagementNM -1.11223  1.21134  0.29730 -0.029662 -0.09172 -0.1446
## ManagementSF -0.19081 -1.03722  0.32552  0.237132  0.15857 -0.0296
## UseHayfield   0.08045  0.76460  0.35627  0.470345  0.07734 -0.1392
## UseHaypastu   0.22445 -0.61214  0.02544  0.013349  0.06101 -0.1916
## UsePasture   -0.48763 -0.02094 -0.53183 -0.668114 -0.20870  0.5131
## Manure1       0.51300  0.02905  0.27896 -0.315483 -0.22424  1.2323
## Manure2       0.72919 -0.03448 -0.70675  0.367543  0.16244 -0.3713
## Manure3      -0.41706 -0.55693 -0.52936  0.221143  0.26656  0.1926
## Manure4       0.44964 -0.96318  1.15039 -0.511103 -0.25361 -0.8552
anova(dune.cca)
#plots
ordiplot(dune.rda, type="text", main="RDA")

ordiplot(dune.cca, type="text", main="CCA")




##Q12

Question 12: Use the results from the CCA ordination in exercise 13 to look at and interpret the variance inflation factor, and to look for warnings about outliers in the explanatory data!

Look at and interpret the inflation factors for each environmental variable! Are there any extreme values?

# VIF for CCA 
dune.cca <- cca(dune ~ . , dune.env)
vif.cca <- vif.cca(dune.cca)
vif.cca
##       A1 Moisture   Manure      Use       BF       HF       NM       SF 
## 1.523245 1.733303 7.878266 1.683073 4.685483 2.945088 7.490125       NA

Function vif.cca gives the variance inflation factors for each constraint or contrast in factor constraints. In partial ordination, conditioning variables are analysed together with constraints. Variance inflation is a diagnostic tool to identify useless constraints. A common rule is that values over 10 indicate redundant constraints.




##Q13

Question 13: Do a CA on the Dune Meadow data set, with the full set of environmental variables included in the analysis, and compare with a CCA on the same data sets.

Compare the constrained and unconstrained ordinations! What is the difference between the methods? When could an unconstrained ordination with explanatory variables be a good option?

# CA with env. variables 
fit <- envfit(dune.ca, dune.env, perm=999)
fit
## 
## ***VECTORS
## 
##               CA1      CA2     r2 Pr(>r)    
## A1        0.99816  0.06061 0.3104  0.046 *  
## Moisture  0.95470  0.29756 0.6802  0.001 ***
## Manure   -0.16779 -0.98582 0.4848  0.003 ** 
## Use       0.73455 -0.67855 0.1616  0.175    
## BF       -0.98156 -0.19113 0.1122  0.311    
## HF       -0.79401 -0.60790 0.1052  0.366    
## NM        0.41219  0.91110 0.7054  0.001 ***
## SF        0.44670 -0.89468 0.2428  0.066 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
ordiplot(dune.ca,  type="text", main="CA with env. variables fitted after")
plot(fit)

#CCA 
dune.cca <- cca(dune ~ . , dune.env.original)
dune.cca
## Call: cca(formula = dune ~ A1 + Moisture + Management + Use + Manure,
## data = dune.env.original)
## 
##               Inertia Proportion Rank
## Total          2.1153     1.0000     
## Constrained    1.5032     0.7106   12
## Unconstrained  0.6121     0.2894    7
## Inertia is scaled Chi-square 
## Some constraints were aliased because they were collinear (redundant)
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3   CCA4   CCA5   CCA6   CCA7   CCA8   CCA9  CCA10  CCA11 
## 0.4671 0.3410 0.1761 0.1532 0.0953 0.0703 0.0589 0.0499 0.0318 0.0260 0.0228 
##  CCA12 
## 0.0108 
## 
## Eigenvalues for unconstrained axes:
##     CA1     CA2     CA3     CA4     CA5     CA6     CA7 
## 0.27237 0.10876 0.08975 0.06305 0.03489 0.02529 0.01798
ordiplot(dune.cca, type="text", main="CA")



#Q14 ### Question 14: Site scores can be calculated in two ways; Weighted Average and Linear Combination. Here you will do both and compare the results. Both options are available in most software, but the default varies.

What is the main difference between the two ways. How to know which method to select? What does the display = “cn” in the last plots mean?

#### Illustration of difference between WA and LC scores

data("dune")
data("dune.env")
dune.env$Moisture<-as.numeric(dune.env$Moisture) #Set Moisture to numeric for the illustration of WA and LC

# Make a CCA with Management and Moisture as explanatory variables
cca.dune<-cca(dune~Management + Moisture, dune.env)
# Look at the result
anova(cca.dune)
summary(eigenvals(cca.dune))
## Importance of components:
##                         CCA1   CCA2    CCA3    CCA4    CA1     CA2     CA3
## Eigenvalue            0.4355 0.2722 0.10304 0.04575 0.3709 0.16307 0.13787
## Proportion Explained  0.2059 0.1287 0.04871 0.02163 0.1754 0.07709 0.06518
## Cumulative Proportion 0.2059 0.3345 0.38326 0.40488 0.5803 0.65734 0.72252
##                           CA4     CA5     CA6     CA7     CA8     CA9    CA10
## Eigenvalue            0.12969 0.09557 0.07909 0.06565 0.05511 0.04546 0.03917
## Proportion Explained  0.06131 0.04518 0.03739 0.03104 0.02605 0.02149 0.01852
## Cumulative Proportion 0.78383 0.82901 0.86640 0.89744 0.92349 0.94498 0.96350
##                          CA11     CA12     CA13     CA14     CA15
## Eigenvalue            0.02712 0.020921 0.016962 0.008008 0.004208
## Proportion Explained  0.01282 0.009891 0.008019 0.003786 0.001989
## Cumulative Proportion 0.97632 0.986206 0.994225 0.998011 1.000000
# Make two plots, one with weighted average (WA) and one with Linear Combination (LC) to calculate the site scores. 
par(mfrow= c(1,2)) # To get two plots beside each other
plot(cca.dune, type = "n", display = c("wa"), scaling = "sites", main = "WA", pch = 3)
points(cca.dune, display = c("wa"), scaling = "sites", pch = 3)
text(cca.dune, display = "wa", scaling = "sites", pos = 2)
ordispider(cca.dune, dune.env$Management, col = "black", label= TRUE, scaling = "sites", display = "wa")

plot(cca.dune, type = "n", display = c("lc"), scaling = "sites", main = "LC", pch = 3)
points(cca.dune, display = c("lc"), scaling = "sites", col = "red", pch = 19)
text(cca.dune, display = "lc", col = "red", scaling = "sites", pos = 3)
ordispider(cca.dune, dune.env$Management, col = "red", label= TRUE, scaling = "sites", display = "lc")

par(mfrow= c(1,1))

# Same as above, but here with sites and species
par(mfrow= c(1,2))
plot(cca.dune, display = c("sp", "wa", "cn"), scaling = "species", main = "Weighted average")
plot(cca.dune, display = c("sp", "lc", "cn"), scaling = "species", main = "Linear combination")

par(mfrow= c(1,1))


##Q15

Question 15: Next, you should investigate whether the four different management regimes (SF standard farming, BF biodynamical farming, HF hobby farming, and NM nature management) have any effect on the species distribution. This is done by using these four management types as environmental variables in a constrained ordination.

Inspect and interpret the results of the permutations and graph! Compare with the CCA on the full set of explanatory data in exercise 13.

To investigate if the vegetational pattern is an effect of management, or if the pattern is random, Monte Carlo permutations will be used.

# CCA using only Management as the environmental variable
dune.ccaMN <- cca(dune ~ Management , dune.env.original)
summary(dune.ccaMN)
## 
## Call:
## cca(formula = dune ~ Management, data = dune.env.original) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.1153     1.0000
## Constrained    0.6038     0.2855
## Unconstrained  1.5114     0.7145
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                         CCA1    CCA2    CCA3    CA1     CA2     CA3     CA4
## Eigenvalue            0.3186 0.18247 0.10274 0.4474 0.20300 0.16301 0.13457
## Proportion Explained  0.1506 0.08626 0.04857 0.2115 0.09597 0.07706 0.06362
## Cumulative Proportion 0.1506 0.23690 0.28547 0.4970 0.59294 0.67000 0.73362
##                           CA5     CA6     CA7     CA8     CA9    CA10   CA11
## Eigenvalue            0.12940 0.09494 0.07904 0.06526 0.05004 0.04321 0.0387
## Proportion Explained  0.06117 0.04488 0.03737 0.03085 0.02366 0.02043 0.0183
## Cumulative Proportion 0.79479 0.83968 0.87704 0.90790 0.93155 0.95198 0.9703
##                          CA12     CA13     CA14     CA15     CA16
## Eigenvalue            0.02385 0.017729 0.009172 0.007959 0.004157
## Proportion Explained  0.01128 0.008382 0.004336 0.003763 0.001965
## Cumulative Proportion 0.98155 0.989936 0.994272 0.998035 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1   CCA2   CCA3
## Eigenvalue            0.3186 0.1825 0.1027
## Proportion Explained  0.5277 0.3022 0.1701
## Cumulative Proportion 0.5277 0.8299 1.0000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##              CCA1     CCA2      CCA3       CA1       CA2      CA3
## Achimill  0.19825  0.72622  0.300280  0.536978 -0.435913  0.64609
## Agrostol -0.10058 -0.73046  0.007614 -0.681142  0.217728  0.03744
## Airaprae -1.87848 -0.05503 -0.068523  2.450771  1.553306  0.96472
## Alopgeni  0.55258 -0.78563  0.076123 -0.233871  0.552230 -0.21144
## Anthodor -0.39419  0.50501 -0.225911  1.160620  0.171184  0.93294
## Bellpere  0.13512  0.17366  0.564850  0.257160 -0.708045 -0.32204
## Bromhord  0.49092  0.60069  0.653072  0.176934 -0.436949  0.13777
## Chenalbu  0.56015 -1.38599  0.350857  0.002875  0.943749 -0.06831
## Cirsarve  0.56015 -1.38599  0.350857  0.405234 -0.657911 -0.47130
## Comapalu -1.87848 -0.05503 -0.068523 -1.746112 -0.512810  0.93566
## Eleopalu -0.70823 -0.37016 -0.119191 -1.553248  0.282390  0.66025
## Elymrepe  0.53962 -0.19035 -0.046524  0.198907 -0.499457 -0.61438
## Empenigr -1.87848 -0.05503 -0.068523  2.309910  2.212605 -0.08448
## Hyporadi -1.36521  0.25217  0.325575  1.904376  1.495023  0.41829
## Juncarti -0.38902  0.03085 -0.512153 -1.205380  0.574275 -0.35684
## Juncbufo  0.55931 -0.45216 -0.375981  0.003634  0.784787 -0.51755
## Lolipere  0.43543  0.19968  0.213581  0.216193 -0.477291 -0.18621
## Planlanc  0.06038  0.66342 -0.325865  0.556636 -0.525653  0.42479
## Poaprat   0.32406  0.11985  0.096677  0.195454 -0.309623 -0.29579
## Poatriv   0.53695 -0.14144 -0.012649  0.046730 -0.006836 -0.07777
## Ranuflam -0.83362 -0.33640 -0.113762 -1.301757  0.264387  0.17231
## Rumeacet  0.55852  0.41250 -1.048979  0.325485 -0.319841  0.40174
## Sagiproc  0.18110 -0.51035  0.108393  0.359930  0.783570 -0.45162
## Salirepe -1.87848 -0.05503 -0.068523  0.123760  0.058476 -0.95975
## Scorautu -0.32939  0.26009  0.120989  0.263833  0.097401 -0.12059
## Trifprat  0.55832  0.63731 -1.223958  0.430250 -0.704454  0.82922
## Trifrepe -0.04953  0.33644  0.186947 -0.090075 -0.056136  0.12455
## Vicilath -0.14617  0.98179  1.261557  0.241368 -0.211654 -0.51063
## Bracruta -0.30221 -0.01389 -0.078786 -0.087945 -0.022522 -0.11640
## Callcusp -1.14689 -0.45432  0.057291 -1.645373 -0.097690  0.92131
## 
## 
## Site scores (weighted averages of species scores)
## 
##        CCA1    CCA2     CCA3       CA1      CA2      CA3
## 1   1.35561  0.4747  1.06564  0.823027 -2.26986 -0.78064
## 2   0.86419  0.8142  1.87245 -0.198170 -0.21094 -0.45453
## 3   0.94050 -0.9209  0.88921  0.342743 -0.72769 -0.44672
## 4   0.77228 -1.1875  1.30265  0.405234 -0.65791 -0.47130
## 5   0.69269  1.5459 -1.41512  0.521056 -0.79020  0.65898
## 6   0.54544  1.7488 -2.38761  0.404007 -0.73194  1.05302
## 7   0.76144  1.5049 -1.05537  0.405051 -0.54999  0.43997
## 8   0.01261 -1.0636 -0.08646 -1.029299  1.29069 -0.66249
## 9   0.76176 -0.7696 -0.94225 -0.400662  0.94009 -1.66619
## 10  0.41252  1.8518  1.64522  0.094693 -0.26410  0.77515
## 11 -0.11176  1.2962  1.56524  0.132854  0.63174 -0.44504
## 12  0.73786 -1.7731 -0.58173  0.207790  1.17765 -0.18459
## 13  0.75284 -1.9690  0.09623  0.002875  0.94375 -0.06831
## 14 -1.85731 -1.0193  0.31713 -1.775732 -0.79413  1.57589
## 15 -1.69911 -1.1108 -0.99568 -1.716492 -0.23149  0.29543
## 16 -1.00521 -2.2417 -0.74574 -1.640220  0.82453  1.87405
## 17 -1.64906  2.1313 -0.06623  2.662063  0.56436  2.53853
## 18 -0.81907  1.1779  0.87653  0.605071 -1.84599 -1.92760
## 19 -2.61902  0.6190  0.42849  2.309910  2.21260 -0.08448
## 20 -2.32115 -1.3309 -1.00044 -1.476716 -0.09132 -0.90420
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##       CCA1     CCA2     CCA3       CA1      CA2      CA3
## 1   0.5601 -1.38599  0.35086  0.823027 -2.26986 -0.78064
## 2   0.4313  1.32739  1.70492 -0.198170 -0.21094 -0.45453
## 3   0.5601 -1.38599  0.35086  0.342743 -0.72769 -0.44672
## 4   0.5601 -1.38599  0.35086  0.405234 -0.65791 -0.47130
## 5   0.5583  0.63731 -1.22396  0.521056 -0.79020  0.65898
## 6   0.5583  0.63731 -1.22396  0.404007 -0.73194  1.05302
## 7   0.5583  0.63731 -1.22396  0.405051 -0.54999  0.43997
## 8   0.5583  0.63731 -1.22396 -1.029299  1.29069 -0.66249
## 9   0.5583  0.63731 -1.22396 -0.400662  0.94009 -1.66619
## 10  0.4313  1.32739  1.70492  0.094693 -0.26410  0.77515
## 11  0.4313  1.32739  1.70492  0.132854  0.63174 -0.44504
## 12  0.5601 -1.38599  0.35086  0.207790  1.17765 -0.18459
## 13  0.5601 -1.38599  0.35086  0.002875  0.94375 -0.06831
## 14 -1.8785 -0.05503 -0.06852 -1.775732 -0.79413  1.57589
## 15 -1.8785 -0.05503 -0.06852 -1.716492 -0.23149  0.29543
## 16  0.5601 -1.38599  0.35086 -1.640220  0.82453  1.87405
## 17 -1.8785 -0.05503 -0.06852  2.662063  0.56436  2.53853
## 18 -1.8785 -0.05503 -0.06852  0.605071 -1.84599 -1.92760
## 19 -1.8785 -0.05503 -0.06852  2.309910  2.21260 -0.08448
## 20 -1.8785 -0.05503 -0.06852 -1.476716 -0.09132 -0.90420
## 
## 
## Biplot scores for constraining variables
## 
##                 CCA1     CCA2     CCA3 CA1 CA2 CA3
## ManagementHF  0.3751  0.42813 -0.82221   0   0   0
## ManagementNM -0.9989 -0.02926 -0.03644   0   0   0
## ManagementSF  0.3648 -0.90262  0.22849   0   0   0
## 
## 
## Centroids for factor constraints
## 
##                 CCA1     CCA2     CCA3 CA1 CA2 CA3
## ManagementBF  0.4313  1.32739  1.70492   0   0   0
## ManagementHF  0.5583  0.63731 -1.22396   0   0   0
## ManagementNM -1.8785 -0.05503 -0.06852   0   0   0
## ManagementSF  0.5601 -1.38599  0.35086   0   0   0
# Permutation test on CCA
MCperm <- permutest(dune.ccaMN, permutations=999)
MCperm
## 
## Permutation test for cca under reduced model 
## 
## Permutation: free
## Number of permutations: 999
##  
## Model: cca(formula = dune ~ Management, data = dune.env.original)
## Permutation test for all constrained eigenvalues
##          Df Inertia      F Pr(>F)   
## Model     3 0.60384 2.1307  0.002 **
## Residual 16 1.51143                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Plot the ordination with only Management
ordiplot(dune.ccaMN, type = "points", main = "CCA only Management")
ordihull(dune.ccaMN,groups = dune.env.original$Management,draw = "polygon",col = "grey70", label = T)




##Q16

Question 16: Next, you want to investigate whether initial soil characteristics (A1 and moisture) also are important for differences in species composition between fields? To do this, do another CCA with A1 and Moisture as environmental variables and perform a significance test.

Explain the results of the permutations! Compare the results in exercises 16 and 17!

# CCA using A1 and moisture as the environmental variable
dune.ccaA1ms <- cca(dune, dune.env[,1:2])
summary(dune.ccaA1ms)
## 
## Call:
## cca(X = dune, Y = dune.env[, 1:2]) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.1153     1.0000
## Constrained    0.5357     0.2533
## Unconstrained  1.5795     0.7467
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## 
## Importance of components:
##                         CCA1    CCA2    CA1    CA2     CA3     CA4     CA5
## Eigenvalue            0.4265 0.10926 0.4044 0.2809 0.17168 0.14665 0.10004
## Proportion Explained  0.2016 0.05165 0.1912 0.1328 0.08116 0.06933 0.04729
## Cumulative Proportion 0.2016 0.25327 0.4444 0.5772 0.65841 0.72774 0.77503
##                           CA6     CA7     CA8     CA9    CA10    CA11    CA12
## Eigenvalue            0.09194 0.07595 0.07326 0.06055 0.04690 0.04554 0.02269
## Proportion Explained  0.04347 0.03590 0.03464 0.02862 0.02217 0.02153 0.01073
## Cumulative Proportion 0.81850 0.85440 0.88904 0.91766 0.93983 0.96136 0.97209
##                           CA13     CA14     CA15     CA16     CA17
## Eigenvalue            0.020748 0.015104 0.010950 0.007971 0.004264
## Proportion Explained  0.009809 0.007141 0.005177 0.003768 0.002016
## Cumulative Proportion 0.981898 0.989039 0.994216 0.997984 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1   CCA2
## Eigenvalue            0.4265 0.1093
## Proportion Explained  0.7961 0.2039
## Cumulative Proportion 0.7961 1.0000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##              CCA1      CCA2      CA1       CA2      CA3      CA4
## Achimill -0.80522 -0.034144  0.38619 -0.039691 -0.22546  0.96607
## Agrostol  0.80865 -0.056511 -0.47652  0.122059  0.14051 -0.06689
## Airaprae  0.45857  0.929764  3.17582 -1.067783  0.42954  0.56399
## Alopgeni  0.48206  0.212415 -0.72292 -0.643671 -0.07254 -0.26939
## Anthodor -0.35975  0.070548  1.37653 -0.370765 -0.43073  0.67331
## Bellpere -0.68384 -0.258107 -0.18712 -0.007862  0.64294  0.10512
## Bromhord -0.73703 -0.006405 -0.18723 -0.075226  0.34288  0.73667
## Chenalbu  1.29869  0.096750 -0.82155 -1.823846 -0.64420 -0.20148
## Cirsarve -0.32362  0.081849 -0.78087 -0.495313  1.57388 -0.13821
## Comapalu  1.87488 -2.431666  0.31732  0.811112  0.65790  0.62535
## Eleopalu  1.40921 -0.388246 -0.16553  1.155390 -0.03457  0.28163
## Elymrepe -0.49671  0.018245 -0.62108 -0.338009  0.53782  0.33711
## Empenigr  0.99750  1.418422  3.18976 -1.141130  0.86374 -0.53939
## Hyporadi  0.23008  0.822734  2.58776 -0.729843  0.67116 -0.15454
## Juncarti  1.06605  0.324950 -0.35744  0.781955 -0.23317 -0.07351
## Juncbufo  0.38686  0.256698 -0.65114 -1.109108 -0.73373 -0.53684
## Lolipere -0.68296  0.121240 -0.23331  0.090110  0.25987  0.11064
## Planlanc -0.88359 -0.341888  0.54345  0.181775 -0.58711 -0.04526
## Poaprat  -0.47409  0.147648 -0.23347 -0.057499  0.22428  0.06748
## Poatriv  -0.12030  0.038928 -0.44519 -0.437623 -0.13930  0.16610
## Ranuflam  1.33049 -0.042806 -0.18729  0.945186 -0.15012  0.05862
## Rumeacet -0.64258 -0.500123  0.06444 -0.221997 -1.25444 -0.10685
## Sagiproc  0.34029  0.378747 -0.02196 -0.817623  0.34366 -0.60087
## Salirepe  0.45049  0.900374  1.06230  1.233516  0.23704 -1.16588
## Scorautu -0.09368 -0.002416  0.42515 -0.010154  0.17886 -0.20059
## Trifprat -0.98927 -0.564091  0.29087  0.151495 -1.58839  0.02149
## Trifrepe -0.03239 -0.209524  0.05308 -0.020468  0.01832  0.13533
## Vicilath -0.90579 -0.038667  0.46915  0.605435  0.77606 -1.05955
## Bracruta  0.05637 -0.129721  0.19033  0.331955 -0.10193 -0.58438
## Callcusp  1.36154 -0.179077 -0.05984  1.635569 -0.03460  0.51697
## 
## 
## Site scores (weighted averages of species scores)
## 
##       CCA1      CCA2      CA1      CA2      CA3      CA4
## 1  -1.2649  0.791185 -0.93676  0.24888  1.48190  1.10604
## 2  -0.8498  0.001097 -0.52316 -0.01984  1.10655  1.04020
## 3  -0.2389  0.419553 -0.88536 -0.49083  0.84021 -0.07468
## 4  -0.1614  0.505509 -0.78087 -0.49531  1.57388 -0.13821
## 5  -1.1319 -1.169583  0.38927 -0.40435 -0.83372  0.44207
## 6  -1.0755 -1.505312  0.37439  0.32710 -1.98034 -0.25573
## 7  -1.0659 -0.672190 -0.01631  0.26832 -1.36315  0.29399
## 8   0.7940  0.411933 -0.62913  0.12266 -0.24099  0.27183
## 9   0.1152  0.609288 -0.86673 -0.69040 -0.44865 -0.12278
## 10 -1.0563 -0.321545  0.22021  0.17724  0.09727  1.45794
## 11 -0.7576  0.451831  0.51559  0.52629  1.08266 -1.56601
## 12  0.5476  0.527784 -0.62514 -1.68048 -0.77125 -1.61783
## 13  0.7085  0.704490 -0.82155 -1.82385 -0.64420 -0.20148
## 14  1.9879 -3.320249  0.25773  1.07992  0.64829  1.44833
## 15  2.0283 -2.735712  0.37691  0.54230  0.66751 -0.19763
## 16  2.0458 -0.760204 -0.59118  1.32873 -0.47799  0.49495
## 17 -0.6410  1.939122  3.15491 -0.95776 -0.22177  2.21905
## 18 -0.5693  0.202225  0.62521  1.19191  0.84167 -2.56413
## 19  0.3779  3.849064  3.18976 -1.14113  0.86374 -0.53939
## 20  1.9395  0.807201  0.04808  2.68327 -0.50175 -0.70282
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##       CCA1     CCA2      CA1      CA2      CA3      CA4
## 1  -1.2002  0.36172 -0.93676  0.24888  1.48190  1.10604
## 2  -1.1086 -0.04053 -0.52316 -0.01984  1.10655  1.04020
## 3  -0.3105  0.02438 -0.88536 -0.49083  0.84021 -0.07468
## 4  -0.3236  0.08185 -0.78087 -0.49531  1.57388 -0.13821
## 5  -0.7419 -1.64952  0.38927 -0.40435 -0.83372  0.44207
## 6  -1.0038 -0.50024  0.37439  0.32710 -1.98034 -0.25573
## 7  -1.2002  0.36172 -0.01631  0.26832 -1.36315  0.29399
## 8   1.0630  1.13110 -0.62913  0.12266 -0.24099  0.27183
## 9   0.3042  0.89380 -0.86673 -0.69040 -0.44865 -0.12278
## 10 -0.4415  0.59902  0.22021  0.17724  0.09727  1.45794
## 11 -1.1086 -0.04053  0.51559  0.52629  1.08266 -1.56601
## 12  0.5792 -0.31295 -0.62514 -1.68048 -0.77125 -1.61783
## 13  1.2987  0.09675 -0.82155 -1.82385 -0.64420 -0.20148
## 14  1.7308 -1.79956  0.25773  1.07992  0.64829  1.44833
## 15  2.0189 -3.06377  0.37691  0.54230  0.66751 -0.19763
## 16  1.2594  0.26914 -0.59118  1.32873 -0.47799  0.49495
## 17 -0.3498  0.19678  3.15491 -0.95776 -0.22177  2.21905
## 18 -0.9645 -0.67263  0.62521  1.19191  0.84167 -2.56413
## 19  0.9975  1.41842  3.18976 -1.14113  0.86374 -0.53939
## 20  0.9713  1.53335  0.04808  2.68327 -0.50175 -0.70282
## 
## 
## Biplot scores for constraining variables
## 
##            CCA1    CCA2 CA1 CA2 CA3 CA4
## A1       0.6034 -0.7974   0   0   0   0
## Moisture 0.9750  0.2222   0   0   0   0
# Permutation test on CCA
MCpermA1ms <- permutest(dune.ccaA1ms, permutations=999)
MCpermA1ms 
## 
## Permutation test for cca under reduced model 
## 
## Permutation: free
## Number of permutations: 999
##  
## Model: cca(X = dune, Y = dune.env[, 1:2])
## Permutation test for all constrained eigenvalues
##          Df Inertia      F Pr(>F)    
## Model     2 0.53573 2.8829  0.001 ***
## Residual 17 1.57954                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Plot the ordination
ordiplot(dune.ccaA1ms, type="text", main="CCA A1 and moisture")




Direct methods with testing of environmental variables

##Q17

Question 17: Partial ordination

Now, you have concluded that both soil characteristics and management type are determining the species composition significantly. But, what if the observed differences between management types are not caused by management type but by initial differences in soil characteristics? To investigate this, you have to test whether there still is a difference in vegetation between management types, after accounting for (i.e. removing) the effects of soil characteristics. This is done by a partial CCA ordination. Partial ordinations are used to eliminate effects of selected variables by specifying them as covariates. In this case, use the two soil characterising variables (A1 and Moisture) as covariates, and the four management types as environmental variables.

What does the permutation test tell you? Also, compare the results with the results from the permutation test on only management types!

#This partials out the effect of A1 and Moisture before analysing the effects of management
MN.cca <- cca(dune ~ Management + Condition(A1 + Moisture), data=dune.env.original)
summary(MN.cca)
## 
## Call:
## cca(formula = dune ~ Management + Condition(A1 + Moisture), data = dune.env.original) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.1153     1.0000
## Conditioned    0.7437     0.3516
## Constrained    0.3954     0.1869
## Unconstrained  0.9761     0.4615
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## after removing the contribution of conditiniong variables
## 
## Importance of components:
##                         CCA1    CCA2    CCA3    CA1     CA2     CA3     CA4
## Eigenvalue            0.2385 0.08653 0.07036 0.3064 0.13191 0.11516 0.10947
## Proportion Explained  0.1739 0.06309 0.05130 0.2234 0.09618 0.08396 0.07982
## Cumulative Proportion 0.1739 0.23701 0.28831 0.5117 0.60787 0.69183 0.77165
##                           CA5     CA6     CA7     CA8     CA9    CA10     CA11
## Eigenvalue            0.07724 0.07575 0.04871 0.03758 0.03106 0.02102 0.012542
## Proportion Explained  0.05632 0.05523 0.03552 0.02740 0.02264 0.01533 0.009145
## Cumulative Proportion 0.82796 0.88320 0.91872 0.94612 0.96876 0.98409 0.993236
##                           CA12
## Eigenvalue            0.009277
## Proportion Explained  0.006764
## Cumulative Proportion 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1    CCA2    CCA3
## Eigenvalue            0.2385 0.08653 0.07036
## Proportion Explained  0.6032 0.21882 0.17794
## Cumulative Proportion 0.6032 0.82206 1.00000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##              CCA1       CCA2      CCA3      CA1        CA2      CA3
## Achimill -0.22184  0.1625650 -0.486532 -0.29895 -0.1910403  0.80559
## Agrostol  0.31801 -0.1591008  0.237934  0.39533 -0.0572108 -0.04906
## Airaprae -1.90497  0.0498869  0.657458 -2.65101 -0.0560655  0.69574
## Alopgeni  0.75017 -0.2904809  0.198070 -0.05494 -0.0179481 -0.14369
## Anthodor -0.67482  0.5664611 -0.023186 -1.13302 -0.4598130  0.53783
## Bellpere -0.11367 -0.2763144 -0.178950  0.22586  0.4032814 -0.05454
## Bromhord  0.11346 -0.0597708 -0.723087  0.06592  0.2219471  0.35933
## Chenalbu  1.65366 -0.8426420  0.255191 -1.21491  0.5419363 -0.08793
## Cirsarve  0.59767 -0.1296937  0.958252  0.37466  0.7316043 -0.84542
## Comapalu -0.99169  0.0882190 -0.518567  0.58438  0.4791343  0.38437
## Eleopalu  0.05742 -0.0080668 -0.116978  0.81272 -0.4722220  0.23229
## Elymrepe  0.29577 -0.0331848  0.246873  0.27977  1.0085693  0.45803
## Empenigr -1.55294 -0.2816661  0.504499 -2.90039  0.7122104 -0.36569
## Hyporadi -1.40452 -0.3342659  0.124456 -2.19156 -0.0002564 -0.14329
## Juncarti -0.11253  0.4647145 -0.179642  0.86164  0.0195413 -0.14119
## Juncbufo  0.43106 -0.1028278  0.177721 -0.30203 -0.0236828  0.02339
## Lolipere  0.21702 -0.0888942 -0.060884  0.16623  0.2453692 -0.06121
## Planlanc -0.28302  0.4634731  0.004689 -0.25092 -0.5717915 -0.09677
## Poaprat   0.13569  0.0008428 -0.009655  0.12030  0.3263968 -0.07370
## Poatriv   0.52883  0.0242477 -0.008592 -0.09088  0.2279744  0.14426
## Ranuflam -0.13726 -0.0519714 -0.010035  0.76566 -0.2142823  0.15268
## Rumeacet  0.26836  0.8823894  0.261872 -0.26470 -0.5038215  0.15287
## Sagiproc  0.22017 -0.2827029  0.171294 -0.54989  0.2153613 -0.62519
## Salirepe -1.69992 -0.4384927  0.770005  0.49821 -0.1961342 -0.48638
## Scorautu -0.39518 -0.0987393 -0.107125 -0.24614  0.0727362 -0.21337
## Trifprat  0.31106  1.1523633  0.368431 -0.26291 -0.8156335  0.12256
## Trifrepe -0.14521 -0.0264189 -0.408072 -0.03594 -0.0125782  0.09669
## Vicilath -0.60359 -0.8681108 -0.984615  0.09323 -0.4824554 -1.41297
## Bracruta -0.26165 -0.0801219  0.147628  0.14077 -0.3547347 -0.48658
## Callcusp -0.41589 -0.3345768  0.110571  1.18584 -0.8150813  1.08169
## 
## 
## Site scores (weighted averages of species scores)
## 
##        CCA1     CCA2    CCA3      CA1      CA2      CA3
## 1   1.15369 -1.27332  0.4310  0.61137  2.18399  1.89365
## 2   0.61229 -1.61462 -1.6552  0.21908  1.29887  1.22585
## 3   0.91451 -0.84774  1.1194  0.39512  0.53538 -0.23449
## 4   0.74161 -0.94528  1.3644  0.37466  0.73160 -0.84542
## 5  -0.08086  1.40948  0.7508 -0.67144  0.18945  0.50159
## 6   0.02924  2.27467  0.8245 -0.24323 -1.32960 -0.14012
## 7   0.39168  1.23604 -0.2425  0.09643 -0.53580  0.40023
## 8   1.25965  0.57148 -0.9681  0.29444  0.79619 -1.07650
## 9  -0.05559  0.67863 -0.4734  0.59314  1.07760  0.29067
## 10 -0.48177  0.90282 -2.7160  0.03464 -0.84210  0.30477
## 11 -0.62527 -2.12839 -0.7924 -0.33410 -0.57320 -2.01847
## 12  0.06671 -0.81435  0.5681 -0.71177 -1.29312 -0.34880
## 13  1.98505 -0.59361  0.3088 -1.21491  0.54194 -0.08793
## 14 -0.61445 -0.66391 -1.6977  0.80540  0.13905  1.83177
## 15 -0.65898  0.29387 -0.4360  0.36335  0.81922 -1.06304
## 16  0.92510  0.06396  0.1286  0.64651 -2.00831  0.86205
## 17 -3.28246  2.50840  0.7073 -2.27695 -1.20848  2.28790
## 18 -1.48400 -2.23671  0.9243  1.00646  0.05868 -1.91972
## 19 -2.94792  0.05320  1.4111 -2.90039  0.71221 -0.36569
## 20 -0.81071  0.06920  1.0103  2.23242 -0.89403  0.30122
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##         CCA1     CCA2    CCA3      CA1      CA2      CA3
## 1   0.799198 -1.62562  1.8039  0.61137  2.18399  1.89365
## 2  -0.004951 -1.34725 -1.5881  0.21908  1.29887  1.22585
## 3   0.606048 -0.12417  0.9430  0.39512  0.53538 -0.23449
## 4   0.597671 -0.12969  0.9583  0.37466  0.73160 -0.84542
## 5   0.469289  1.25664  0.0800 -0.67144  0.18945  0.50159
## 6   0.301753  1.14623  0.3854 -0.24323 -1.32960 -0.14012
## 7   0.176100  1.06342  0.6144  0.09643 -0.53580  0.40023
## 8   0.879780  1.74702 -0.6594  0.29444  0.79619 -1.07650
## 9  -0.363187  1.16959 -0.3949  0.59314  1.07760  0.29067
## 10 -0.340508  0.06035 -2.1894  0.03464 -0.84210  0.30477
## 11 -0.004951 -1.34725 -1.5881 -0.33410 -0.57320 -2.01847
## 12  0.435824 -1.40351  0.4739 -0.71177 -1.29312 -0.34880
## 13  1.653661 -0.84264  0.2552 -1.21491  0.54194 -0.08793
## 14 -1.083840  0.02749 -0.3506  0.80540  0.13905  1.83177
## 15 -0.899549  0.14895 -0.6865  0.36335  0.81922 -1.06304
## 16  1.628531 -0.85920  0.3010  0.64651 -2.00831  0.86205
## 17 -2.433019  0.54722  0.8869 -2.27695 -1.20848  2.28790
## 18 -2.063955 -0.83829  1.4271  1.00646  0.05868 -1.91972
## 19 -1.552942 -0.28167  0.5045 -2.90039  0.71221 -0.36569
## 20 -1.569696 -0.29271  0.5350  2.23242 -0.89403  0.30122
## 
## 
## Biplot scores for constraining variables
## 
##                 CCA1     CCA2    CCA3 CA1 CA2 CA3
## ManagementHF  0.1944  0.85341 0.01122   0   0   0
## ManagementNM -0.8301 -0.09912 0.21079   0   0   0
## ManagementSF  0.6036 -0.46400 0.47327   0   0   0
## 
## 
## Centroids for factor constraints
## 
##                 CCA1    CCA2    CCA3 CA1 CA2 CA3
## ManagementBF -0.1283 -0.8299 -1.8091   0   0   0
## ManagementHF  0.2894  1.2704  0.0167   0   0   0
## ManagementNM -1.5611 -0.1864  0.3964   0   0   0
## ManagementSF  0.9269 -0.7125  0.7267   0   0   0
#permutation test
permutest(MN.cca, permutations=999)
## 
## Permutation test for cca under reduced model 
## 
## Permutation: free
## Number of permutations: 999
##  
## Model: cca(formula = dune ~ Management + Condition(A1 + Moisture), data
## = dune.env.original)
## Permutation test for all constrained eigenvalues
##          Df Inertia      F Pr(>F)  
## Model     3 0.39543 1.6205  0.021 *
## Residual 12 0.97610                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Plot the ordination
ordiplot(MN.cca, type = "points")
ordihull(MN.cca,groups = dune.env.original$Management,draw = "polygon",col = "grey70", label = T)




##Q18

Question 18: The amount of applied manure is a third important variable given in the environmental data set associated with the dune meadow data set. Repeat exercise 17, but with manure added to the set of covariates.

Inspect and interpret the results of the permutations!

#This partials out the effect of A1, Moisture, and Manure before analysing the effects of management
MN.cca2 <- cca(dune ~ Management + Condition(A1 + Moisture + Manure), data=dune.env.original)
summary(MN.cca2)
## 
## Call:
## cca(formula = dune ~ Management + Condition(A1 + Moisture + Manure),      data = dune.env.original) 
## 
## Partitioning of scaled Chi-square:
##               Inertia Proportion
## Total          2.1153    1.00000
## Conditioned    1.2207    0.57708
## Constrained    0.1524    0.07204
## Unconstrained  0.7422    0.35088
## 
## Eigenvalues, and their contribution to the scaled Chi-square 
## after removing the contribution of conditiniong variables
## 
## Importance of components:
##                         CCA1    CCA2    CA1    CA2    CA3     CA4     CA5
## Eigenvalue            0.1019 0.05047 0.2759 0.1130 0.1047 0.07949 0.06436
## Proportion Explained  0.1139 0.05642 0.3084 0.1263 0.1170 0.08885 0.07195
## Cumulative Proportion 0.1139 0.17035 0.4787 0.6051 0.7221 0.81095 0.88290
##                           CA6     CA7     CA8     CA9
## Eigenvalue            0.04147 0.03105 0.01842 0.01381
## Proportion Explained  0.04636 0.03471 0.02059 0.01543
## Cumulative Proportion 0.92926 0.96397 0.98457 1.00000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         CCA1    CCA2
## Eigenvalue            0.1019 0.05047
## Proportion Explained  0.6688 0.33119
## Cumulative Proportion 0.6688 1.00000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## 
## 
## Species scores
## 
##               CCA1      CCA2      CA1        CA2       CA3      CA4
## Achimill -0.133335  0.120821 -0.11336  0.6510292 -0.337080  0.22393
## Agrostol -0.059233 -0.132637  0.37444 -0.0370997  0.095375 -0.07786
## Airaprae  0.407915 -0.230470 -2.44958  1.0332387  0.418358  0.36362
## Alopgeni -0.280823 -0.253038 -0.01362 -0.2317652  0.041553  0.10891
## Anthodor  0.639819 -0.006334 -0.83017  0.5633922 -0.000667  0.16112
## Bellpere -0.434577  0.069176  0.21427 -0.3756168 -0.309418  0.31025
## Bromhord -0.599725  0.321928  0.10185  0.0261899 -0.332718 -0.04073
## Chenalbu  0.057009 -1.809598 -0.99123 -1.3005988 -2.653184  0.02770
## Cirsarve  0.055897  0.173517  0.05198 -0.8024478  0.546950 -0.82929
## Comapalu -0.596782  0.942411  0.16931  0.5922871 -0.074486 -0.77810
## Eleopalu -0.128649  0.070611  0.82317  0.5764876  0.621879  0.16780
## Elymrepe -0.006479  0.312140 -0.04221  0.3913889 -0.497683  0.46659
## Empenigr  0.356332 -0.139879 -2.96848 -0.1566239  0.522892 -0.37091
## Hyporadi  0.131937 -0.381471 -2.18339  0.4041744  0.743463 -0.02112
## Juncarti  0.168285  0.541592  0.75315  0.0785039  0.371408  0.73033
## Juncbufo -0.059236 -0.346550 -0.24256 -0.2636189 -0.567659 -0.18092
## Lolipere -0.124276  0.092667 -0.01431  0.0985747  0.046798 -0.02100
## Planlanc  0.515751 -0.077257 -0.06547 -0.0155724  0.129069 -0.06399
## Poaprat  -0.080901  0.116929 -0.04072 -0.0256861 -0.125095  0.11124
## Poatriv  -0.056841 -0.011211 -0.08430 -0.1019155 -0.424954  0.07988
## Ranuflam -0.039091  0.005367  0.81474  0.0492154 -0.156240  0.03827
## Rumeacet  0.974169  0.193476 -0.00166  0.0550963  0.027334 -0.19856
## Sagiproc -0.127410 -0.181867 -0.63145 -0.5495718  0.316077 -0.35896
## Salirepe  0.248228 -0.408801  0.66350 -0.7497857 -0.070381  0.23967
## Scorautu -0.144092  0.019065 -0.26341 -0.2057033  0.063209 -0.02494
## Trifprat  1.300548  0.240744  0.09341  0.0640604  0.167503 -0.48798
## Trifrepe -0.190514  0.114486 -0.02230  0.0004746 -0.178668 -0.46336
## Vicilath -0.479706 -0.825859 -0.01391 -0.8227623  0.713571 -0.16779
## Bracruta  0.143858 -0.215431  0.20953 -0.3155810  0.462924  0.10961
## Callcusp -0.030780 -0.303185  1.35836  1.1001438 -0.074428 -1.09279
## 
## 
## Site scores (weighted averages of species scores)
## 
##         CCA1    CCA2      CA1      CA2      CA3      CA4
## 1   0.076849  1.6054 -0.38254  2.52074 -1.15972  1.05976
## 2  -2.949104  1.3252  0.12543  0.38912 -0.76331  0.55160
## 3   0.000495 -0.5491  0.11367 -0.23158 -0.09345  0.45606
## 4  -0.031180 -0.1541  0.05198 -0.80245  0.54695 -0.82929
## 5   1.275980  0.1262 -0.37104  0.08451 -0.62045  0.96726
## 6   1.992850  0.1218  0.33239 -0.07570  0.55582 -0.86651
## 7   0.458162  0.9887 -0.03962  0.39302 -0.01533 -0.99690
## 8  -0.341512  1.4118 -0.09208 -0.80160  0.81681  0.41771
## 9   0.634800  1.5936  0.12543  0.38912 -0.76331  0.55160
## 10 -0.403563  0.6018  0.42285  0.07216 -0.57670 -0.06757
## 11 -0.290887 -2.9002 -0.73283 -0.60769  1.77679 -0.63318
## 12 -0.761760 -1.9123 -0.15051 -0.46695  0.91598 -0.66193
## 13 -0.306982 -2.4614 -0.99123 -1.30060 -2.65318  0.02770
## 14 -0.984376  1.0251  0.54638  1.60362 -1.15108 -3.50187
## 15  0.113636  1.6486 -0.20777 -0.41904  1.00211  1.94566
## 16  0.165589 -0.4483  1.15086  1.79585  1.68170  0.67435
## 17  1.249098  0.2014 -1.67122  2.81803  0.26156  1.46543
## 18 -1.372702 -1.5766  0.98715 -2.14783 -0.12260  0.66276
## 19  1.211626 -1.3634 -2.96848 -0.15662  0.52289 -0.37091
## 20  0.057337  0.6223  2.64850 -0.26686 -0.39501  0.35216
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##         CCA1     CCA2      CA1      CA2      CA3      CA4
## 1  -0.232348 -0.85528 -0.38254  2.52074 -1.15972  1.05976
## 2  -2.714892  0.72966  0.12543  0.38912 -0.76331  0.55160
## 3   0.041672  0.18967  0.11367 -0.23158 -0.09345  0.45606
## 4   0.055897  0.17352  0.05198 -0.80245  0.54695 -0.82929
## 5   1.504037  0.29428 -0.37104  0.08451 -0.62045  0.96726
## 6   1.788549 -0.02879  0.33239 -0.07570  0.55582 -0.86651
## 7  -0.122942  0.86104 -0.03962  0.39302 -0.01533 -0.99690
## 8  -0.006331  2.16478 -0.09208 -0.80160  0.81681  0.41771
## 9   0.869012  0.99804  0.12543  0.38912 -0.76331  0.55160
## 10 -0.266549 -0.23023  0.42285  0.07216 -0.57670 -0.06757
## 11 -0.782403 -1.00056 -0.73283 -0.60769  1.77679 -0.63318
## 12 -1.042814 -1.19765 -0.15051 -0.46695  0.91598 -0.66193
## 13  0.057009 -1.80960 -0.99123 -1.30060 -2.65318  0.02770
## 14 -0.440301  0.76472  0.54638  1.60362 -1.15108 -3.50187
## 15 -0.753264  1.12010 -0.20777 -0.41904  1.00211  1.94566
## 16  0.099686 -1.85806  1.15086  1.79585  1.68170  0.67435
## 17  0.485289 -0.36636 -1.67122  2.81803  0.26156  1.46543
## 18 -0.087468 -1.07208  0.98715 -2.14783 -0.12260  0.66276
## 19  0.356332 -0.13988 -2.96848 -0.15662  0.52289 -0.37091
## 20  0.384783 -0.17219  2.64850 -0.26686 -0.39501  0.35216
## 
## 
## Biplot scores for constraining variables
## 
##                 CCA1    CCA2 CA1 CA2 CA3 CA4
## ManagementHF  0.5735  0.5495   0   0   0   0
## ManagementSF -0.1000 -0.5202   0   0   0   0
## 
## 
## Centroids for factor constraints
## 
##                 CCA1     CCA2 CA1 CA2 CA3 CA4
## ManagementBF -1.2865 -0.09634   0   0   0   0
## ManagementHF  0.8538  0.81795   0   0   0   0
## ManagementSF -0.1536 -0.79878   0   0   0   0
#permutation test
permutest(MN.cca2, permutations=999)
## 
## Permutation test for cca under reduced model 
## 
## Permutation: free
## Number of permutations: 999
##  
## Model: cca(formula = dune ~ Management + Condition(A1 + Moisture +
## Manure), data = dune.env.original)
## Permutation test for all constrained eigenvalues
##          Df Inertia      F Pr(>F)
## Model     2 0.15239 0.9239  0.542
## Residual  9 0.74220




##Q19

Question 19: Forward selection

We have now concluded that management type is of subordinate importance. Now, you want to test the importance or significance of individual environmental variables. This is done by Forward selection. What is the effect of Forward selection? In what way do these results differ from the results obtained when the whole Environmental data table was used as explanatory variables? Were the set of explanatory variables selected by the Forward selection the best set? Why would one perform a p-value correction?

# A first CCA with all environmental variables
cca.dune.all <- cca(dune ~ ., data = dune.env) # the dot after the tilde (~) means "include all from data"
anova(cca.dune.all) # Test overall significance of the explanatory variables. If not significant no need to proceed with step wise selection
adjR2.cca.dune <- RsquareAdj(cca.dune.all)$adj.r.squared # How much the adjusted R2 explains when all variables in dune.env are included
adjR2.cca.dune
## [1] 0.2299496
# A null model with only intercept, needed for the step wise selection
cca.dune.0 <- cca(dune ~ 1, data = dune.env) # model containing only species matrix and intercept

# The step wise selection, based on R2
sel.env.dune <- ordiR2step(cca.dune.0, scope = formula (cca.dune.all), R2scope = adjR2.cca.dune, direction = 'forward', permutations = 999)
## Step: R2.adj= 0 
## Call: dune ~ 1 
##  
##                 R2.adjusted
## <All variables>  0.23030613
## + Management     0.15184053
## + Moisture       0.15084705
## + Manure         0.10001447
## + A1             0.05470899
## + Use            0.01489235
## <none>           0.00000000
## 
##              Df    AIC      F Pr(>F)   
## + Management  3 86.935 2.1307  0.004 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.1523214 
## Call: dune ~ Management 
##  
##                 R2.adjusted
## + Moisture        0.2462733
## <All variables>   0.2303061
## + A1              0.2000052
## + Use             0.1532069
## <none>            0.1523214
## + Manure          0.1466277
sel.env.dune
## Call: cca(formula = dune ~ Management, data = dune.env)
## 
##               Inertia Proportion Rank
## Total          2.1153     1.0000     
## Constrained    0.6038     0.2855    3
## Unconstrained  1.5114     0.7145   16
## Inertia is scaled Chi-square 
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3 
## 0.3186 0.1825 0.1027 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8    CA9   CA10   CA11 
## 0.4474 0.2030 0.1630 0.1346 0.1294 0.0949 0.0790 0.0653 0.0500 0.0432 0.0387 
##   CA12   CA13   CA14   CA15   CA16 
## 0.0239 0.0177 0.0092 0.0080 0.0042
sel.env.dune$anova

Since there is (a potentially high) number of tests of significance during the forward selection procedure, it is better to apply a correction for multiple testing issue. Here we apply the Holm correction.

sel.env.dune_adj <- sel.env.dune
sel.env.dune_adj$anova$`Pr(>F)` <- p.adjust(sel.env.dune$anova$`Pr(>F)`, method = 'holm', n = ncol(dune.env)) # 7 other methods available
sel.env.dune_adj$anova




##Q20

Question 20: Interpreting ordination results using species traits

Experienced plant ecologists may already have looked at the species in the Dune Meadow ordination graphs and concluded that species with similar traits occur together. This is possible if you have good knowledge about plant species ecological preferences, and if there are relatively few species in your dataset. In this exercise we will use tabulated data on species ecological preferences (the Ellenberg indicator values) to interpret the results of the ordinations. Which Ellenberg values are most important for the distribution of the species? What do the different axes represent in terms of environmental gradients?

#load Ellenberg data
dune.ell <- readRDS("Ellenberg.RDS")
dune.mean.ell <- readRDS("Mean_Ellenberg.RDS")
#The vegan implementation of forward selection of variables does not accept NA values.
#To get around this we will replace the two NAs with the mean value for all plots in that column.
#Do not do this in a real analysis without seriously thinking about how it affects your results!
## Code to replace missing values with the column mean, and save the results as a new dataframe
## called dune.mean.ell.impute (since we are imputing the missing values).
dune.mean.ell.impute <- data.frame(
    sapply(
        dune.mean.ell,
        function(x) ifelse(is.na(x),
            mean(x, na.rm = TRUE),
            x)))
#CCA analysis
#create global model with CCA (including all variables) and test it's significance, and if it is significant,
#we use the ordistep function with appropriate arguments to do forward selection of variables.
cca1 <- cca(dune ~ ., data = dune.mean.ell.impute) # full model (with all explanatory variables)
anova(cca1) #overall model is significant
# Start by making ordinations
cca0 <- cca (dune~ 1, data = dune.mean.ell.impute) # empty model only with intercept
cca1 <- cca(dune ~ ., data = dune.mean.ell.impute, na.action = na.omit) # full model (with all explanatory variables)
cca1
## Call: cca(formula = dune ~ L + T + K + F + R + N + S, data =
## dune.mean.ell.impute, na.action = na.omit)
## 
##               Inertia Proportion Rank
## Total          2.1153     1.0000     
## Constrained    1.3313     0.6294    7
## Unconstrained  0.7840     0.3706   12
## Inertia is scaled Chi-square 
## 
## Eigenvalues for constrained axes:
##   CCA1   CCA2   CCA3   CCA4   CCA5   CCA6   CCA7 
## 0.4646 0.3298 0.1774 0.1490 0.1091 0.0597 0.0417 
## 
## Eigenvalues for unconstrained axes:
##     CA1     CA2     CA3     CA4     CA5     CA6     CA7     CA8     CA9    CA10 
## 0.24495 0.14477 0.08921 0.07974 0.06167 0.04827 0.03879 0.02701 0.01952 0.01271 
##    CA11    CA12 
## 0.01006 0.00729
plot(cca1, main="CCA, all data")

#Stepwise approach, using "ordistep"
step2<-ordistep(cca0, scope = formula(cca1), direction="forward")
## 
## Start: dune ~ 1 
## 
##     Df    AIC      F Pr(>F)   
## + N  1 85.551 4.1029  0.005 **
## + L  1 86.138 3.4630  0.005 **
## + F  1 86.635 2.9367  0.005 **
## + T  1 86.624 2.9479  0.010 **
## + S  1 87.111 2.4440  0.015 * 
## + R  1 87.551 1.9994  0.025 * 
## + K  1 87.850 1.7028  0.080 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dune ~ N 
## 
##     Df    AIC      F Pr(>F)   
## + F  1 83.748 3.5595  0.005 **
## + S  1 84.490 2.8114  0.005 **
## + L  1 84.501 2.8000  0.005 **
## + R  1 84.753 2.5526  0.015 * 
## + T  1 85.234 2.0874  0.020 * 
## + K  1 85.992 1.3778  0.160   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dune ~ N + F 
## 
##     Df    AIC      F Pr(>F)  
## + T  1 82.828 2.5154  0.020 *
## + R  1 83.824 1.6159  0.060 .
## + S  1 83.715 1.7119  0.075 .
## + L  1 83.543 1.8649  0.080 .
## + K  1 83.881 1.5660  0.100 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dune ~ N + F + T 
## 
##     Df    AIC      F Pr(>F)  
## + L  1 82.381 1.9525  0.015 *
## + S  1 82.728 1.6611  0.065 .
## + R  1 82.582 1.7830  0.070 .
## + K  1 83.052 1.3932  0.170  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dune ~ N + F + T + L 
## 
##     Df    AIC      F Pr(>F)  
## + S  1 81.835 1.9003  0.055 .
## + R  1 82.299 1.5355  0.120  
## + K  1 82.328 1.5135  0.145  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
step2$anova # which three are most significant in the final selection?




##Q21

Question 21: Decomposition of variance

Perform the analyses and interpret the results! How much of the total variation (% of All) is explained by: - uniquely by Management (effect of Soil removed), - uniquely by Soil (effect of Management removed), - jointly by Soil and Management (the interaction) - Not by Soil or Management?

In usual analysis of variance experiments, the variance is decomposed into components. The same can be done in multivariate analyses. In this exercise you will decompose the variance in the Dune meadow data set into different variance components. You will use two groups of variance components: Group 1 is Management and Group 2 is the soil variables (A1 and Moisture). In ordinations, the variance is expressed by the sum of the eigenvalues.

The result gives you the fraction of the total variation that is explained by: a) uniquely by Management (effect of soil removed), b) uniquely by soil (effect of Management removed), c) jointly by Soil and Management (the interaction).

The total explained variation is the sum of a), b) and c).

The total variation in the dataset is the sum of all unconstrained eigenvalues.

data("dune")
data("dune.env")
#as we did in Question 3:PCA, we need to reformulate the environmental data to make
#new columns for each of the management methods.
dummy_management <- as.data.frame(model.matrix( ~ Management - 1, data=dune.env )) 
#add these to the dataset
dune.env <- dune.env %>% select(A1, Moisture, Manure, Use) %>% cbind(.,dummy_management) 
dune.env$Moisture <- as.numeric(as.character(dune.env$Moisture)) #make numeric
dune.env$Manure <- as.numeric(as.character(dune.env$Manure))
dune.env$Use <- as.numeric(dune.env$Use)
#make column names shorter
dune.env <- dune.env %>% rename(BF = ManagementBF, HF = ManagementHF,
                                NM = ManagementNM, SF = ManagementSF)
## variance partitioning
management <- dune.env[,c("BF", "HF", "NM", "SF")]
soil <- dune.env[,c("A1", "Moisture")]
# examine the explanatory variable of each class of variables.
varp <- varpart(dune, management, soil)
varp
## 
## Partition of variance in RDA 
## 
## Call: varpart(Y = dune, X = management, soil)
## 
## Explanatory tables:
## X1:  management
## X2:  soil 
## 
## No. of explanatory tables: 2 
## Total variation (SS): 1598.4 
##             Variance: 84.124 
## No. of observations: 20 
## 
## Partition table:
##                      Df R.squared Adj.R.squared Testable
## [a+b] = X1            3   0.34747       0.22512     TRUE
## [b+c] = X2            2   0.25771       0.17038     TRUE
## [a+b+c] = X1+X2       5   0.51552       0.34249     TRUE
## Individual fractions                                    
## [a] = X1|X2           3                 0.17211     TRUE
## [b]                   0                 0.05301    FALSE
## [c] = X2|X1           2                 0.11737     TRUE
## [d] = Residuals                         0.65751    FALSE
## ---
## Use function 'rda' to test significance of fractions of interest
plot (varp, digits = 2, Xnames = c('Management', 'Soil'), bg = c('red', 'blue'))




##Q22

Question 22: Analysing a time series with vegetation data

Is there a significant change in species composition over time? What is the effect of defining Plots as cofactors?

Quite often vegetation ecologists have data from repeated inventories in permanent plots. The overall question to answer is if there has been a significant and directional change in the species composition. In this exercise we have rearranged the Dune Meadow data so that it consists of only 10 plots, each analysed 2 times (same species as in all other exercises, just grouping and rearrangement of plots). We want to analyse if there has been a consistent change in species composition between the two inventories. We are thus not interested in differences between the ten plots.

#This partials out the effect of Plot before analysing the effects of Time
#First load the time series versions of the data
SpeTS <- readRDS("SpeTS.RDS")
EnvTS <- readRDS("EnvTS.RDS")
time.cca <- cca(SpeTS ~ Time + Condition(Plot), data=EnvTS)
time.cca
## Call: cca(formula = SpeTS ~ Time + Condition(Plot), data = EnvTS)
## 
##               Inertia Proportion Rank
## Total         2.11526    1.00000     
## Conditional   0.20908    0.09885    1
## Constrained   0.15941    0.07536    1
## Unconstrained 1.74677    0.82579   17
## Inertia is scaled Chi-square 
## 
## Eigenvalues for constrained axes:
##    CCA1 
## 0.15941 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.4215 0.2761 0.2461 0.1722 0.1358 0.0999 0.0913 0.0783 
## (Showing 8 of 17 unconstrained eigenvalues)
treat <- EnvTS$Time #
colvec <- c("red2", "green4") #set colours to be applied to different levels of factor "Time"
plot(time.cca, type = "n", display = c("lc"))#plots the axes
with(time.cca, points(time.cca, display = c("lc"), col = colvec[treat],#plots the points
                      pch = 21, xlim = c(-2,2), bg = colvec[treat]))

#permutation test
with(EnvTS, anova(time.cca, by="term", perm=500, strata=Plot))
## Set of permutations < 'minperm'. Generating entire set.
#Compare with results for time as only explanatory variable and no cofactors
time2.cca <- cca(SpeTS ~ Time, data=EnvTS)
time2.cca
## Call: cca(formula = SpeTS ~ Time, data = EnvTS)
## 
##               Inertia Proportion Rank
## Total         2.11526    1.00000     
## Constrained   0.14240    0.06732    1
## Unconstrained 1.97287    0.93268   18
## Inertia is scaled Chi-square 
## 
## Eigenvalues for constrained axes:
##   CCA1 
## 0.1424 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8 
## 0.4796 0.3721 0.2530 0.1736 0.1365 0.1063 0.0921 0.0783 
## (Showing 8 of 18 unconstrained eigenvalues)
#permutation test
with(EnvTS, anova(time2.cca, by="term", perm=500))
## plot ellipsoid hulls
treat <- EnvTS$Time
plot(time2.cca, type = "n")#plots the axes
with(time2.cca, points(time2.cca, col = colvec[treat],#plots the points
                      pch = 21, bg = colvec[treat]))
ordihull(time2.cca,groups=treat,draw="polygon",col="grey70",label=T)#draw hulls around area of each level of factor "Time"




##Q23

Question 23: A multivariate Before-After-Control- Impact (BACI) study

Did the treatment have an effect? What was permuted?

In this exercise we will continue to use the rearranged Dune Meadow data from exercise 24. A difference is that the plots are now divided into four groups: 1. Control plots before a treatment (i.e. not treated) 2. Control plots after a treatment (i.e. not treated) 3. Impact plots before treatment (i.e. not treated) 4. Impact plots after treatment (this is where the treatment is made)

The four groups are indicated in variable Treat in the environmental data set (EnvTS). The question you want to answer is if the treatment caused a change in species composition that is significantly different from the change in the control plots.

baci.cca <- cca(SpeTS ~ Treat + Condition(Plot + Time), data=EnvTS)
baci.cca
## Call: cca(formula = SpeTS ~ Treat + Condition(Plot + Time), data =
## EnvTS)
## 
##               Inertia Proportion Rank
## Total          2.1153     1.0000     
## Conditional    0.3685     0.1742    2
## Constrained    0.2443     0.1155    1
## Unconstrained  1.5024     0.7103   16
## Inertia is scaled Chi-square 
## 
## Eigenvalues for constrained axes:
##    CCA1 
## 0.24432 
## 
## Eigenvalues for unconstrained axes:
##    CA1    CA2    CA3    CA4    CA5    CA6    CA7    CA8    CA9   CA10   CA11 
## 0.3925 0.2472 0.1740 0.1453 0.1182 0.0933 0.0829 0.0570 0.0532 0.0441 0.0324 
##   CA12   CA13   CA14   CA15   CA16 
## 0.0218 0.0141 0.0115 0.0082 0.0069
with(EnvTS, anova(baci.cca, by="term", perm=500, strata=Treat))
#Here with() is a special function that makes variables in dune.env visible to
#the following command. If you only type Moisture in an R prompt, you will get
#an error of missing variables
## plot results
treat <- EnvTS$Treat
ordiplot(baci.cca,type="points")
ordihull(baci.cca,groups=treat,draw="polygon",col="grey70",label=T)#draw hulls around area of each level of factor "Time"




##Q24

Question 24: Procustes rotation and Co-correspondence analysis

Are the two data sets (species and bugs) correlated? Give a very short ecological interpretation of the results!

These are two ways of comparing ordinations or the structure in related multivariate datasets. They are used to compare how objects differ when they are described by different sets of descriptors, or if there are repeated samplings of the same objects using the same descriptors.

library(cocorresp)
# DCA on the plant species data 
DCASpecies <- decorana(dune)
# DCA on the insect data (called Bugs)
Bugs <- readRDS("Bugs.RDS")
DCABugs <- decorana(Bugs)
DCABugs
## 
## Call:
## decorana(veg = Bugs) 
## 
## Detrended correspondence analysis with 26 segments.
## Rescaling of axes with 4 iterations.
## 
##                   DCA1   DCA2    DCA3   DCA4
## Eigenvalues     0.6785 0.3236 0.17269 0.1292
## Decorana values 0.6853 0.2470 0.03452 0.0207
## Axis lengths    4.2552 2.2877 1.63208 1.6331
# procrustes analysis on the two DCAs from different datasets
ProC <- procrustes(DCASpecies, DCABugs)

plot(ProC, kind=1, type="text")

# Labels show the position of the samples in the second ordination, and arrows point to their positions in the target ordination
# Function protest tests the non-randomness (`significance') between two configurations
ProCsig <- protest(dune.dca, DCABugs)
ProCsig
## 
## Call:
## protest(X = dune.dca, Y = DCABugs) 
## 
## Procrustes Sum of Squares (m12 squared):        0.4492 
## Correlation in a symmetric Procrustes rotation: 0.7422 
## Significance:  0.001 
## 
## Permutation: free
## Number of permutations: 999
# Co-correspondence analysis
CoCor <- coca(dune ~ ., data = Bugs, method = "symmetric")
summary(CoCor)
## 
## Symmetric Co-Correspondence Analysis
## 
## Call: symcoca(y = y, x = x, n.axes = n.axes, R0 = weights, symmetric =
## symmetric, nam.dat = nam.dat)
## 
## Inertia:
##       Total Explained Residual
## dune: 2.115     1.672    0.444
## Bugs: 2.358     2.280    0.078
## 
## Eigenvalues:
##    COCA 1     COCA 2     COCA 3     COCA 4     COCA 5     COCA 6     COCA 7  
## 0.2751016  0.2551218  0.0662534  0.0210179  0.0075231  0.0035365  0.0014639  
##    COCA 8     COCA 9  
## 0.0008311  0.0002375
corAxis(CoCor)
##    COCA 1    COCA 2    COCA 3    COCA 4    COCA 5    COCA 6    COCA 7    COCA 8 
## 0.9021743 0.9103017 0.8259081 0.8133787 0.8165470 0.7043982 0.7022924 0.5065464 
##    COCA 9 
## 0.7012526




##Q25

Question 25 : Classification

How do the different clustering techniques differ?

Start by testing the non-hierarchical K-means clustering, using Multivar / K-Means. Note that it is recommended that K-means need at least 100 observations (500 according to some sources) to be reliable! Let us ignore the sample size issue for the moment, and ask for 4 clusters (for later comparison with the four management types). Test different hierarchical agglomeration algorithms and similarity indices. Use at least the Euclidian and the Bray-Curtis similarity measures for the hierarchical clustering technique. Try also the Ward’s method!

library(dendextend)
#make a distance matrix using the Bray-Curtis method
dune.dist <- vegdist(dune, method = "bray")
# K-Means Cluster Analysis
fit <- kmeans(dune, 4, nstart = 50) # we specify a 4 cluster solution
#As the final result of k-means clustering result is sensitive to the random starting assignments, we specify nstart = 50. This means #that R will try 50 different random starting assignments and then select the best results corresponding to the one with the lowest #within cluster variation.
#Note that k-means clustering involves randomness so you won't neccessarily get exactly the same results if you repeat it.
#There are also several different algorithms that can be used in k-means clustering, and different software 
#have can use different defaults.Check out the help file for the kmeans function if you want to know more about these.
#cluster number for each data point
fit$cluster
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 
##  3  3  2  2  3  3  3  2  2  3  1  2  2  4  4  4  1  1  1  4
#cluster sizes
fit$size
## [1] 4 6 6 4
#compare clusters to management category
table(fit$cluster, dune.env.original$Management)
##    
##     BF HF NM SF
##   1  1  0  3  0
##   2  0  2  0  4
##   3  2  3  0  1
##   4  0  0  3  1
#try some approaches to hierarchical clustering
dend2 <- dune %>% # data
        dist(method = "euclidean") %>% # calculate a distance matrix, choose method 
        hclust(method = "ward.D") %>% # Hierarchical clustering, choose method 
        as.dendrogram # Turn the object into a dendrogram.
dend2 %>% set("labels_color") %>% set("branches_lwd") %>%  set("branches_k_color", k = 4) %>% plot

#the dist() function does not include Bray-Curtis but we have already made a distance matrix using 
#the vegdist() function from vegan.
dend3 <- dune.dist %>% # data (Bray-Curtis)
        hclust(method = "aver") %>% # Hierarchical clustering, choose method 
        as.dendrogram # Turn the object into a dendrogram.
dend3 %>% set("labels_color") %>% set("branches_lwd") %>%  set("branches_k_color", k = 4) %>% plot

#we can colour the labels based on their grouping in the kmeans analysis to compare clusters
dend3 %>% set("labels_color", fit$cluster) %>% set("branches_lwd") %>%  set("branches_k_color", k = 4) %>% plot

#tanglegram is a nice function to compare dendrograms that you might want to try out as an extra
tanglegram(dend2, dend3)



#Q26 ### Question 26: Classification 2

Classification 2

#to follow

##Q27

Question 27: Principal Response Curves

PRC

#to follow

##Q28

Question 28: ANOSIM & SIMPER

What does the result tell you? What other types of predefined clusters could you use?

ANOSIM (ANalysis Of Similarities) is a non-parametric test of significant difference between two or more groups, based on any distance measure. In this case, use the clusters from the K-means exercise. Change to Bray-Curtis similarity index.

dune.dist <- vegdist(dune, method = "bray") #create distance matrix based on Bray-Curtis method
dune.ano <- anosim(dune.dist, fit$cluster) #compare groupings based on cluster analysis
summary(dune.ano)
## 
## Call:
## anosim(x = dune.dist, grouping = fit$cluster) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.8226 
##       Significance: 0.001 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.116 0.171 0.208 0.263 
## 
## Dissimilarity ranks between and within classes:
##           0%    25%    50%    75%  100%   N
## Between 12.5 76.875 115.50 153.50 188.0 148
## 1       10.0 62.250  66.25 101.75 134.0   6
## 2        3.0 14.250  22.00  29.50  53.5  15
## 3        1.0  6.000  23.50  55.75  91.0  15
## 4        6.0 15.250  19.50  32.00  60.0   6
#plot(dune.ano)
#title(sub = "Anosim, groups = cluster analysis")
#ANOSIM gives you the P value and a R value. R value close to 1 indicates high separation between levels of your factor while R value #close to 0 indicate no separation between levels of your factor.
dune.ano2 <- anosim(dune.dist, dune.env.original$Management) #try comparing groupings based on management instead
summary(dune.ano2)
## 
## Call:
## anosim(x = dune.dist, grouping = dune.env.original$Management) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.2579 
##       Significance: 0.008 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.122 0.153 0.198 0.230 
## 
## Dissimilarity ranks between and within classes:
##         0%   25%    50%     75%  100%   N
## Between  4 58.50 104.00 145.500 188.0 147
## BF       5 15.25  25.50  41.250  57.0   3
## HF       1  7.25  46.25  68.125  89.5  10
## NM       6 64.75 124.50 156.250 181.0  15
## SF       3 32.75  53.50  99.250 184.0  15
#plot(dune.ano2, sub = "Anosim, groups = management")

SIMPER

Compare the different groups and check which species that contributes mostly to the difference between pairs of clusters.

After a significant ANOSIM, you may want to know which species are primarily responsible for the observed difference between clusters. SIMPER (Similarity Percentage) will do this for you. The test does not come with significance testing. In the output tables, taxa are sorted in descending order of contribution to group difference.

sim <- simper(dune, fit$cluster) #start with kmeans clusters
summary(sim)
## 
## Contrast: 3_2 
## 
##           average       sd  ratio    ava    avb cumsum
## Alopgeni 0.062264 0.033663 1.8496 0.3333 5.0000 0.1062
## Agrostol 0.060675 0.021567 2.8133 0.0000 4.6667 0.2097
## Lolipere 0.042051 0.035413 1.1875 5.3333 2.8333 0.2814
## Planlanc 0.036403 0.027638 1.3171 3.0000 0.0000 0.3435
## Sagiproc 0.032522 0.021454 1.5159 0.0000 2.5000 0.3990
## Elymrepe 0.030536 0.028642 1.0662 2.0000 2.3333 0.4511
## Achimill 0.029361 0.010739 2.7340 2.3333 0.0000 0.5012
## Rumeacet 0.028649 0.025240 1.1351 2.3333 0.6667 0.5501
## Trifrepe 0.027416 0.017993 1.5237 3.3333 2.1667 0.5968
## Poatriv  0.026724 0.027452 0.9735 4.6667 5.5000 0.6424
## Anthodor 0.026238 0.020621 1.2724 2.1667 0.0000 0.6872
## Juncbufo 0.024559 0.024039 1.0216 0.3333 1.8333 0.7291
## Bromhord 0.024328 0.019807 1.2283 2.0000 0.5000 0.7706
## Bracruta 0.021988 0.021354 1.0297 2.0000 2.0000 0.8081
## Scorautu 0.018939 0.013653 1.3872 2.8333 2.1667 0.8404
## Trifprat 0.017870 0.021259 0.8406 1.5000 0.0000 0.8709
## Poaprat  0.017839 0.019132 0.9324 3.5000 3.1667 0.9013
## Juncarti 0.016978 0.024765 0.6856 0.0000 1.3333 0.9303
## Bellpere 0.014809 0.014581 1.0156 1.1667 0.6667 0.9555
## Ranuflam 0.009039 0.013245 0.6824 0.0000 0.6667 0.9709
## Eleopalu 0.008600 0.019773 0.4349 0.0000 0.6667 0.9856
## Cirsarve 0.004034 0.009258 0.4357 0.0000 0.3333 0.9925
## Chenalbu 0.002369 0.005468 0.4333 0.0000 0.1667 0.9965
## Vicilath 0.002033 0.004618 0.4403 0.1667 0.0000 1.0000
## Airaprae 0.000000 0.000000    NaN 0.0000 0.0000 1.0000
## Comapalu 0.000000 0.000000    NaN 0.0000 0.0000 1.0000
## Empenigr 0.000000 0.000000    NaN 0.0000 0.0000 1.0000
## Hyporadi 0.000000 0.000000    NaN 0.0000 0.0000 1.0000
## Salirepe 0.000000 0.000000    NaN 0.0000 0.0000 1.0000
## Callcusp 0.000000 0.000000    NaN 0.0000 0.0000 1.0000
## 
## Contrast: 3_1 
## 
##           average      sd  ratio    ava  avb cumsum
## Poatriv  0.071134 0.02222 3.2013 4.6667 0.00 0.1104
## Lolipere 0.065842 0.04953 1.3292 5.3333 2.25 0.2126
## Bracruta 0.042351 0.03280 1.2912 2.0000 3.25 0.2783
## Scorautu 0.037921 0.03176 1.1939 2.8333 4.50 0.3372
## Planlanc 0.036828 0.02270 1.6226 3.0000 2.00 0.3944
## Trifrepe 0.036651 0.02824 1.2979 3.3333 1.75 0.4512
## Hyporadi 0.035220 0.02812 1.2524 0.0000 2.25 0.5059
## Elymrepe 0.035093 0.03850 0.9114 2.0000 0.00 0.5604
## Rumeacet 0.033387 0.03621 0.9221 2.3333 0.00 0.6122
## Anthodor 0.032373 0.03200 1.0118 2.1667 2.00 0.6624
## Poaprat  0.030417 0.02631 1.1560 3.5000 2.00 0.7097
## Bromhord 0.029558 0.02489 1.1873 2.0000 0.00 0.7555
## Achimill 0.028319 0.01536 1.8441 2.3333 0.50 0.7995
## Salirepe 0.022678 0.02398 0.9458 0.0000 1.50 0.8347
## Trifprat 0.021280 0.02555 0.8330 1.5000 0.00 0.8677
## Airaprae 0.020683 0.02216 0.9333 0.0000 1.25 0.8998
## Sagiproc 0.018199 0.01995 0.9123 0.0000 1.25 0.9281
## Bellpere 0.017753 0.01816 0.9775 1.1667 0.50 0.9556
## Vicilath 0.011242 0.01211 0.9287 0.1667 0.75 0.9731
## Empenigr 0.007323 0.01323 0.5536 0.0000 0.50 0.9844
## Juncbufo 0.005090 0.01172 0.4345 0.3333 0.00 0.9923
## Alopgeni 0.004937 0.01136 0.4347 0.3333 0.00 1.0000
## Agrostol 0.000000 0.00000    NaN 0.0000 0.00 1.0000
## Chenalbu 0.000000 0.00000    NaN 0.0000 0.00 1.0000
## Cirsarve 0.000000 0.00000    NaN 0.0000 0.00 1.0000
## Comapalu 0.000000 0.00000    NaN 0.0000 0.00 1.0000
## Eleopalu 0.000000 0.00000    NaN 0.0000 0.00 1.0000
## Juncarti 0.000000 0.00000    NaN 0.0000 0.00 1.0000
## Ranuflam 0.000000 0.00000    NaN 0.0000 0.00 1.0000
## Callcusp 0.000000 0.00000    NaN 0.0000 0.00 1.0000
## 
## Contrast: 3_4 
## 
##           average       sd  ratio    ava  avb  cumsum
## Lolipere 0.084704 0.038558 2.1968 5.3333 0.00 0.09603
## Eleopalu 0.080308 0.026995 2.9749 0.0000 5.25 0.18707
## Agrostol 0.076247 0.020401 3.7374 0.0000 5.00 0.27351
## Poatriv  0.061858 0.026316 2.3506 4.6667 0.50 0.34364
## Poaprat  0.054987 0.019512 2.8181 3.5000 0.00 0.40598
## Trifrepe 0.044678 0.034197 1.3065 3.3333 1.75 0.45663
## Planlanc 0.042303 0.032421 1.3048 3.0000 0.00 0.50458
## Bracruta 0.038932 0.025172 1.5466 2.0000 3.00 0.54872
## Ranuflam 0.038285 0.014364 2.6654 0.0000 2.50 0.59212
## Callcusp 0.038035 0.025820 1.4731 0.0000 2.50 0.63524
## Juncarti 0.037813 0.024529 1.5416 0.0000 2.50 0.67811
## Achimill 0.034414 0.012390 2.7776 2.3333 0.00 0.71712
## Elymrepe 0.033755 0.036292 0.9301 2.0000 0.00 0.75539
## Rumeacet 0.032477 0.034996 0.9280 2.3333 0.00 0.79221
## Anthodor 0.030478 0.024163 1.2614 2.1667 0.00 0.82676
## Bromhord 0.028722 0.023955 1.1990 2.0000 0.00 0.85932
## Scorautu 0.027794 0.017188 1.6170 2.8333 1.50 0.89083
## Trifprat 0.020706 0.024733 0.8372 1.5000 0.00 0.91431
## Salirepe 0.018308 0.033071 0.5536 0.0000 1.25 0.93506
## Alopgeni 0.016792 0.024399 0.6882 0.3333 1.00 0.95410
## Bellpere 0.016654 0.017807 0.9352 1.1667 0.00 0.97298
## Comapalu 0.016530 0.017597 0.9393 0.0000 1.00 0.99172
## Juncbufo 0.004940 0.011311 0.4367 0.3333 0.00 0.99732
## Vicilath 0.002365 0.005413 0.4368 0.1667 0.00 1.00000
## Airaprae 0.000000 0.000000    NaN 0.0000 0.00 1.00000
## Chenalbu 0.000000 0.000000    NaN 0.0000 0.00 1.00000
## Cirsarve 0.000000 0.000000    NaN 0.0000 0.00 1.00000
## Empenigr 0.000000 0.000000    NaN 0.0000 0.00 1.00000
## Hyporadi 0.000000 0.000000    NaN 0.0000 0.00 1.00000
## Sagiproc 0.000000 0.000000    NaN 0.0000 0.00 1.00000
## 
## Contrast: 2_1 
## 
##           average       sd  ratio    ava  avb cumsum
## Poatriv  0.086332 0.034365 2.5122 5.5000 0.00 0.1160
## Alopgeni 0.079013 0.037490 2.1076 5.0000 0.00 0.2222
## Agrostol 0.071974 0.024236 2.9697 4.6667 0.00 0.3189
## Lolipere 0.045515 0.035550 1.2803 2.8333 2.25 0.3801
## Bracruta 0.036395 0.024788 1.4683 2.0000 3.25 0.4290
## Scorautu 0.035223 0.020888 1.6863 2.1667 4.50 0.4763
## Elymrepe 0.034442 0.037075 0.9290 2.3333 0.00 0.5226
## Hyporadi 0.034208 0.026214 1.3050 0.0000 2.25 0.5686
## Anthodor 0.032868 0.034303 0.9582 0.0000 2.00 0.6127
## Poaprat  0.032448 0.023611 1.3743 3.1667 2.00 0.6563
## Planlanc 0.031236 0.018827 1.6591 0.0000 2.00 0.6983
## Sagiproc 0.030263 0.024672 1.2266 2.5000 1.25 0.7390
## Juncbufo 0.029473 0.030854 0.9552 1.8333 0.00 0.7786
## Salirepe 0.022104 0.022682 0.9745 0.0000 1.50 0.8083
## Juncarti 0.020055 0.029245 0.6858 1.3333 0.00 0.8352
## Airaprae 0.020009 0.020637 0.9696 0.0000 1.25 0.8621
## Trifrepe 0.016639 0.017363 0.9583 2.1667 1.75 0.8845
## Bellpere 0.012598 0.015378 0.8192 0.6667 0.50 0.9014
## Vicilath 0.010842 0.012071 0.8982 0.0000 0.75 0.9160
## Ranuflam 0.010799 0.015811 0.6830 0.6667 0.00 0.9305
## Rumeacet 0.010455 0.015293 0.6836 0.6667 0.00 0.9445
## Eleopalu 0.010180 0.023430 0.4345 0.6667 0.00 0.9582
## Achimill 0.009284 0.016490 0.5630 0.0000 0.50 0.9707
## Empenigr 0.007150 0.012680 0.5639 0.0000 0.50 0.9803
## Bromhord 0.007088 0.016294 0.4350 0.5000 0.00 0.9898
## Cirsarve 0.004725 0.010863 0.4350 0.3333 0.00 0.9962
## Chenalbu 0.002855 0.006584 0.4336 0.1667 0.00 1.0000
## Comapalu 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## Trifprat 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## Callcusp 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## 
## Contrast: 2_4 
## 
##           average       sd  ratio    ava  avb cumsum
## Poatriv  0.076730 0.035297 2.1738 5.5000 0.50 0.1165
## Eleopalu 0.068316 0.032100 2.1282 0.6667 5.25 0.2202
## Alopgeni 0.065917 0.039379 1.6739 5.0000 1.00 0.3203
## Poaprat  0.046536 0.024779 1.8781 3.1667 0.00 0.3910
## Lolipere 0.040993 0.034867 1.1757 2.8333 0.00 0.4532
## Sagiproc 0.037422 0.023982 1.5605 2.5000 0.00 0.5101
## Callcusp 0.037049 0.023751 1.5599 0.0000 2.50 0.5663
## Trifrepe 0.035339 0.018937 1.8661 2.1667 1.75 0.6200
## Elymrepe 0.033468 0.035748 0.9362 2.3333 0.00 0.6708
## Juncarti 0.032485 0.024290 1.3373 1.3333 2.50 0.7201
## Bracruta 0.030206 0.019581 1.5426 2.0000 3.00 0.7660
## Juncbufo 0.028542 0.029595 0.9644 1.8333 0.00 0.8093
## Ranuflam 0.026823 0.018649 1.4383 0.6667 2.50 0.8501
## Agrostol 0.021487 0.019803 1.0851 4.6667 5.00 0.8827
## Salirepe 0.017875 0.031700 0.5639 0.0000 1.25 0.9098
## Comapalu 0.016027 0.016444 0.9746 0.0000 1.00 0.9342
## Rumeacet 0.010135 0.014719 0.6885 0.6667 0.00 0.9496
## Bellpere 0.009538 0.013832 0.6896 0.6667 0.00 0.9641
## Scorautu 0.009421 0.013442 0.7008 2.1667 1.50 0.9784
## Bromhord 0.006897 0.015787 0.4369 0.5000 0.00 0.9888
## Cirsarve 0.004598 0.010525 0.4369 0.3333 0.00 0.9958
## Chenalbu 0.002757 0.006317 0.4365 0.1667 0.00 1.0000
## Achimill 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## Airaprae 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## Anthodor 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## Empenigr 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## Hyporadi 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## Planlanc 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## Trifprat 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## Vicilath 0.000000 0.000000    NaN 0.0000 0.00 1.0000
## 
## Contrast: 1_4 
## 
##           average      sd  ratio  ava  avb cumsum
## Eleopalu 0.098265 0.03006 3.2686 0.00 5.25 0.1177
## Agrostol 0.093231 0.02119 4.4001 0.00 5.00 0.2294
## Scorautu 0.052080 0.02919 1.7844 4.50 1.50 0.2917
## Ranuflam 0.046847 0.01632 2.8700 0.00 2.50 0.3479
## Callcusp 0.046477 0.03138 1.4810 0.00 2.50 0.4035
## Juncarti 0.046157 0.02950 1.5647 0.00 2.50 0.4588
## Trifrepe 0.042804 0.03664 1.1683 1.75 1.75 0.5101
## Bracruta 0.042329 0.03965 1.0675 3.25 3.00 0.5608
## Hyporadi 0.041620 0.03195 1.3028 2.25 0.00 0.6106
## Anthodor 0.040746 0.04350 0.9366 2.00 0.00 0.6594
## Lolipere 0.038632 0.04988 0.7744 2.25 0.00 0.7057
## Planlanc 0.038219 0.02329 1.6412 2.00 0.00 0.7515
## Poaprat  0.036516 0.02753 1.3265 2.00 0.00 0.7952
## Salirepe 0.036288 0.03376 1.0749 1.50 1.25 0.8387
## Airaprae 0.024651 0.02570 0.9594 1.25 0.00 0.8682
## Sagiproc 0.021247 0.02299 0.9241 1.25 0.00 0.8936
## Comapalu 0.020538 0.02173 0.9453 0.00 1.00 0.9182
## Alopgeni 0.017127 0.03097 0.5530 0.00 1.00 0.9388
## Vicilath 0.013007 0.01462 0.8898 0.75 0.00 0.9543
## Achimill 0.011816 0.02128 0.5553 0.50 0.00 0.9685
## Bellpere 0.009189 0.01651 0.5567 0.50 0.00 0.9795
## Poatriv  0.008564 0.01549 0.5530 0.00 0.50 0.9898
## Empenigr 0.008557 0.01536 0.5570 0.50 0.00 1.0000
## Bromhord 0.000000 0.00000    NaN 0.00 0.00 1.0000
## Chenalbu 0.000000 0.00000    NaN 0.00 0.00 1.0000
## Cirsarve 0.000000 0.00000    NaN 0.00 0.00 1.0000
## Elymrepe 0.000000 0.00000    NaN 0.00 0.00 1.0000
## Juncbufo 0.000000 0.00000    NaN 0.00 0.00 1.0000
## Rumeacet 0.000000 0.00000    NaN 0.00 0.00 1.0000
## Trifprat 0.000000 0.00000    NaN 0.00 0.00 1.0000
## Permutation: free
## Number of permutations: 0
sim2 <- simper(dune, dune.env.original$Management) #try management groupings
summary(sim2)
## 
## Contrast: SF_BF 
## 
##           average       sd  ratio    ava    avb  cumsum
## Agrostol 0.061374 0.034193 1.7949 4.6667 0.0000 0.09824
## Alopgeni 0.052667 0.036476 1.4439 4.3333 0.6667 0.18255
## Lolipere 0.048116 0.039445 1.2198 3.0000 6.0000 0.25957
## Trifrepe 0.046297 0.025525 1.8138 1.3333 4.6667 0.33368
## Poatriv  0.046020 0.033801 1.3615 4.6667 3.6667 0.40734
## Scorautu 0.043697 0.024922 1.7534 1.3333 4.3333 0.47729
## Bromhord 0.033677 0.025860 1.3023 0.5000 2.6667 0.53120
## Achimill 0.030152 0.020821 1.4482 0.1667 2.3333 0.57947
## Planlanc 0.028585 0.021549 1.3265 0.0000 2.0000 0.62522
## Elymrepe 0.028074 0.029778 0.9428 2.0000 1.3333 0.67016
## Bracruta 0.025501 0.023902 1.0669 2.0000 2.0000 0.71098
## Poaprat  0.025129 0.023967 1.0485 2.5000 4.0000 0.75121
## Sagiproc 0.024326 0.022149 1.0983 1.8333 0.6667 0.79014
## Bellpere 0.019859 0.017088 1.1622 0.6667 1.6667 0.82193
## Eleopalu 0.018611 0.042958 0.4333 1.3333 0.0000 0.85172
## Anthodor 0.017543 0.025804 0.6798 0.0000 1.3333 0.87981
## Juncbufo 0.016031 0.023708 0.6762 1.1667 0.0000 0.90547
## Vicilath 0.014671 0.013306 1.1026 0.0000 1.0000 0.92895
## Hyporadi 0.010286 0.015198 0.6768 0.0000 0.6667 0.94542
## Ranuflam 0.009306 0.013595 0.6845 0.6667 0.0000 0.96031
## Juncarti 0.006979 0.016109 0.4333 0.5000 0.0000 0.97148
## Callcusp 0.006979 0.016109 0.4333 0.5000 0.0000 0.98266
## Rumeacet 0.004526 0.010444 0.4333 0.3333 0.0000 0.98990
## Cirsarve 0.003983 0.009185 0.4336 0.3333 0.0000 0.99628
## Chenalbu 0.002326 0.005370 0.4333 0.1667 0.0000 1.00000
## Airaprae 0.000000 0.000000    NaN 0.0000 0.0000 1.00000
## Comapalu 0.000000 0.000000    NaN 0.0000 0.0000 1.00000
## Empenigr 0.000000 0.000000    NaN 0.0000 0.0000 1.00000
## Salirepe 0.000000 0.000000    NaN 0.0000 0.0000 1.00000
## Trifprat 0.000000 0.000000    NaN 0.0000 0.0000 1.00000
## 
## Contrast: SF_HF 
## 
##           average       sd  ratio    ava avb  cumsum
## Agrostol 0.047380 0.031273 1.5151 4.6667 1.4 0.08351
## Alopgeni 0.046433 0.032897 1.4115 4.3333 1.6 0.16535
## Lolipere 0.041986 0.027007 1.5546 3.0000 4.0 0.23935
## Planlanc 0.039198 0.033208 1.1804 0.0000 3.0 0.30844
## Rumeacet 0.038992 0.027369 1.4247 0.3333 3.2 0.37716
## Elymrepe 0.031877 0.029550 1.0787 2.0000 2.0 0.43334
## Poatriv  0.028466 0.021522 1.3227 4.6667 4.8 0.48352
## Bracruta 0.025261 0.021044 1.2004 2.0000 2.8 0.52804
## Eleopalu 0.024974 0.038877 0.6424 1.3333 0.8 0.57206
## Poaprat  0.023932 0.019180 1.2478 2.5000 3.4 0.61424
## Anthodor 0.023409 0.021430 1.0923 0.0000 1.8 0.65550
## Sagiproc 0.023144 0.020479 1.1301 1.8333 0.8 0.69629
## Trifprat 0.023080 0.023432 0.9850 0.0000 1.8 0.73697
## Juncarti 0.022850 0.025677 0.8899 0.5000 1.6 0.77724
## Trifrepe 0.022383 0.019487 1.1486 1.3333 2.8 0.81669
## Juncbufo 0.021643 0.022237 0.9733 1.1667 1.2 0.85484
## Scorautu 0.020509 0.016422 1.2489 1.3333 2.8 0.89099
## Achimill 0.015183 0.011393 1.3326 0.1667 1.2 0.91775
## Bromhord 0.013375 0.014504 0.9222 0.5000 0.8 0.94132
## Ranuflam 0.010661 0.013387 0.7964 0.6667 0.4 0.96011
## Bellpere 0.009991 0.012571 0.7948 0.6667 0.4 0.97772
## Callcusp 0.006623 0.015076 0.4393 0.5000 0.0 0.98939
## Cirsarve 0.003809 0.008669 0.4394 0.3333 0.0 0.99611
## Chenalbu 0.002208 0.005025 0.4393 0.1667 0.0 1.00000
## Airaprae 0.000000 0.000000    NaN 0.0000 0.0 1.00000
## Comapalu 0.000000 0.000000    NaN 0.0000 0.0 1.00000
## Empenigr 0.000000 0.000000    NaN 0.0000 0.0 1.00000
## Hyporadi 0.000000 0.000000    NaN 0.0000 0.0 1.00000
## Salirepe 0.000000 0.000000    NaN 0.0000 0.0 1.00000
## Vicilath 0.000000 0.000000    NaN 0.0000 0.0 1.00000
## 
## Contrast: SF_NM 
## 
##           average       sd  ratio    ava    avb cumsum
## Poatriv  0.078284 0.040947 1.9118 4.6667 0.0000 0.1014
## Alopgeni 0.071219 0.046958 1.5167 4.3333 0.0000 0.1936
## Agrostol 0.056508 0.044176 1.2792 4.6667 2.1667 0.2667
## Lolipere 0.054851 0.059914 0.9155 3.0000 0.3333 0.3378
## Eleopalu 0.048027 0.047168 1.0182 1.3333 2.1667 0.3999
## Poaprat  0.040724 0.031790 1.2810 2.5000 0.6667 0.4527
## Bracruta 0.040008 0.034398 1.1631 2.0000 2.8333 0.5045
## Elymrepe 0.035598 0.038515 0.9243 2.0000 0.0000 0.5506
## Scorautu 0.032475 0.034813 0.9328 1.3333 3.1667 0.5926
## Trifrepe 0.030430 0.031634 0.9619 1.3333 1.8333 0.6320
## Sagiproc 0.030304 0.030477 0.9943 1.8333 0.5000 0.6712
## Salirepe 0.029275 0.032014 0.9144 0.0000 1.8333 0.7092
## Anthodor 0.024541 0.036694 0.6688 0.0000 1.3333 0.7409
## Callcusp 0.022763 0.029443 0.7731 0.5000 1.1667 0.7704
## Ranuflam 0.022566 0.022819 0.9889 0.6667 1.3333 0.7996
## Juncarti 0.022543 0.028598 0.7883 0.5000 1.1667 0.8288
## Hyporadi 0.020108 0.031291 0.6426 0.0000 1.1667 0.8548
## Juncbufo 0.019860 0.029034 0.6840 1.1667 0.0000 0.8806
## Planlanc 0.015420 0.022769 0.6772 0.0000 0.8333 0.9005
## Airaprae 0.014883 0.021881 0.6802 0.0000 0.8333 0.9198
## Bellpere 0.012317 0.015921 0.7737 0.6667 0.3333 0.9357
## Comapalu 0.011883 0.017407 0.6826 0.0000 0.6667 0.9511
## Achimill 0.009294 0.014931 0.6224 0.1667 0.3333 0.9632
## Bromhord 0.007172 0.016333 0.4391 0.5000 0.0000 0.9724
## Rumeacet 0.005590 0.012751 0.4384 0.3333 0.0000 0.9797
## Empenigr 0.005225 0.012001 0.4354 0.0000 0.3333 0.9864
## Cirsarve 0.004782 0.010889 0.4391 0.3333 0.0000 0.9926
## Chenalbu 0.002893 0.006602 0.4382 0.1667 0.0000 0.9964
## Vicilath 0.002792 0.006425 0.4345 0.0000 0.1667 1.0000
## Trifprat 0.000000 0.000000    NaN 0.0000 0.0000 1.0000
## 
## Contrast: BF_HF 
## 
##           average      sd  ratio    ava avb  cumsum
## Rumeacet 0.038666 0.02606 1.4838 0.0000 3.2 0.08163
## Poatriv  0.033301 0.02579 1.2911 3.6667 4.8 0.15194
## Planlanc 0.031852 0.01830 1.7401 2.0000 3.0 0.21918
## Bromhord 0.028651 0.01799 1.5926 2.6667 0.8 0.27967
## Lolipere 0.028431 0.02215 1.2834 6.0000 4.0 0.33970
## Elymrepe 0.027822 0.02959 0.9404 1.3333 2.0 0.39843
## Trifrepe 0.025838 0.01656 1.5603 4.6667 2.8 0.45298
## Anthodor 0.023582 0.02042 1.1547 1.3333 1.8 0.50277
## Achimill 0.023426 0.01474 1.5893 2.3333 1.2 0.55223
## Bracruta 0.022733 0.01802 1.2617 2.0000 2.8 0.60022
## Alopgeni 0.021610 0.02308 0.9363 0.6667 1.6 0.64584
## Trifprat 0.021514 0.02207 0.9747 0.0000 1.8 0.69126
## Juncarti 0.020084 0.02555 0.7860 0.0000 1.6 0.73367
## Scorautu 0.019318 0.01356 1.4241 4.3333 2.8 0.77445
## Bellpere 0.018290 0.01486 1.2305 1.6667 0.4 0.81306
## Agrostol 0.017605 0.02284 0.7708 0.0000 1.4 0.85023
## Juncbufo 0.015000 0.02066 0.7260 0.0000 1.2 0.88190
## Vicilath 0.012848 0.01140 1.1274 1.0000 0.0 0.90903
## Sagiproc 0.011685 0.01297 0.9008 0.6667 0.8 0.93369
## Eleopalu 0.010169 0.02111 0.4817 0.0000 0.8 0.95516
## Hyporadi 0.008950 0.01312 0.6824 0.6667 0.0 0.97406
## Poaprat  0.007203 0.01010 0.7133 4.0000 3.4 0.98927
## Ranuflam 0.005084 0.01055 0.4817 0.0000 0.4 1.00000
## Airaprae 0.000000 0.00000    NaN 0.0000 0.0 1.00000
## Chenalbu 0.000000 0.00000    NaN 0.0000 0.0 1.00000
## Cirsarve 0.000000 0.00000    NaN 0.0000 0.0 1.00000
## Comapalu 0.000000 0.00000    NaN 0.0000 0.0 1.00000
## Empenigr 0.000000 0.00000    NaN 0.0000 0.0 1.00000
## Salirepe 0.000000 0.00000    NaN 0.0000 0.0 1.00000
## Callcusp 0.000000 0.00000    NaN 0.0000 0.0 1.00000
## 
## Contrast: BF_NM 
## 
##           average      sd  ratio    ava    avb cumsum
## Lolipere 0.090681 0.02644 3.4292 6.0000 0.3333 0.1243
## Poatriv  0.054684 0.04465 1.2248 3.6667 0.0000 0.1992
## Poaprat  0.052511 0.01813 2.8966 4.0000 0.6667 0.2712
## Trifrepe 0.051287 0.02756 1.8611 4.6667 1.8333 0.3415
## Bromhord 0.039689 0.02920 1.3590 2.6667 0.0000 0.3959
## Bracruta 0.035723 0.02869 1.2452 2.0000 2.8333 0.4448
## Eleopalu 0.033759 0.03573 0.9449 0.0000 2.1667 0.4911
## Agrostol 0.033446 0.03473 0.9630 0.0000 2.1667 0.5369
## Achimill 0.033190 0.02338 1.4198 2.3333 0.3333 0.5824
## Scorautu 0.031356 0.02026 1.5480 4.3333 3.1667 0.6254
## Anthodor 0.028060 0.03295 0.8517 1.3333 1.3333 0.6638
## Planlanc 0.027319 0.02193 1.2458 2.0000 0.8333 0.7013
## Salirepe 0.026770 0.02927 0.9145 0.0000 1.8333 0.7379
## Bellpere 0.023529 0.01909 1.2322 1.6667 0.3333 0.7702
## Hyporadi 0.021721 0.02450 0.8864 0.6667 1.1667 0.8000
## Ranuflam 0.020314 0.02275 0.8928 0.0000 1.3333 0.8278
## Elymrepe 0.019993 0.02926 0.6833 1.3333 0.0000 0.8552
## Callcusp 0.017833 0.02681 0.6653 0.0000 1.1667 0.8796
## Juncarti 0.017694 0.02600 0.6806 0.0000 1.1667 0.9039
## Vicilath 0.015773 0.01447 1.0902 1.0000 0.1667 0.9255
## Sagiproc 0.015432 0.01857 0.8310 0.6667 0.5000 0.9466
## Airaprae 0.013410 0.01969 0.6809 0.0000 0.8333 0.9650
## Comapalu 0.010739 0.01571 0.6835 0.0000 0.6667 0.9797
## Alopgeni 0.009997 0.01463 0.6833 0.6667 0.0000 0.9934
## Empenigr 0.004787 0.01105 0.4332 0.0000 0.3333 1.0000
## Chenalbu 0.000000 0.00000    NaN 0.0000 0.0000 1.0000
## Cirsarve 0.000000 0.00000    NaN 0.0000 0.0000 1.0000
## Juncbufo 0.000000 0.00000    NaN 0.0000 0.0000 1.0000
## Rumeacet 0.000000 0.00000    NaN 0.0000 0.0000 1.0000
## Trifprat 0.000000 0.00000    NaN 0.0000 0.0000 1.0000
## 
## Contrast: HF_NM 
## 
##           average       sd  ratio ava    avb  cumsum
## Poatriv  0.071553 0.013681 5.2302 4.8 0.0000 0.09913
## Lolipere 0.054533 0.029625 1.8408 4.0 0.3333 0.17468
## Rumeacet 0.046546 0.030806 1.5109 3.2 0.0000 0.23917
## Poaprat  0.041750 0.018852 2.2146 3.4 0.6667 0.29701
## Planlanc 0.041633 0.029560 1.4084 3.0 0.8333 0.35469
## Bracruta 0.035340 0.020104 1.7579 2.8 2.8333 0.40365
## Eleopalu 0.032043 0.032315 0.9916 0.8 2.1667 0.44805
## Agrostol 0.031915 0.028887 1.1048 1.4 2.1667 0.49227
## Trifrepe 0.030372 0.022871 1.3280 2.8 1.8333 0.53434
## Elymrepe 0.029811 0.038676 0.7708 2.0 0.0000 0.57565
## Anthodor 0.028717 0.024799 1.1580 1.8 1.3333 0.61543
## Juncarti 0.027414 0.028537 0.9607 1.6 1.1667 0.65341
## Trifprat 0.025843 0.025972 0.9951 1.8 0.0000 0.68922
## Salirepe 0.025340 0.027291 0.9285 0.0 1.8333 0.72432
## Alopgeni 0.024459 0.032399 0.7549 1.6 0.0000 0.75821
## Scorautu 0.020703 0.014125 1.4658 2.8 3.1667 0.78689
## Ranuflam 0.019285 0.019939 0.9672 0.4 1.3333 0.81361
## Juncbufo 0.018181 0.024648 0.7376 1.2 0.0000 0.83880
## Hyporadi 0.017141 0.026548 0.6457 0.0 1.1667 0.86255
## Callcusp 0.016833 0.024901 0.6760 0.0 1.1667 0.88587
## Achimill 0.016555 0.014900 1.1111 1.2 0.3333 0.90881
## Sagiproc 0.015282 0.016535 0.9243 0.8 0.5000 0.92998
## Airaprae 0.012605 0.018243 0.6910 0.0 0.8333 0.94744
## Bromhord 0.012094 0.015169 0.7973 0.8 0.0000 0.96420
## Comapalu 0.010105 0.014556 0.6942 0.0 0.6667 0.97820
## Bellpere 0.008801 0.013732 0.6409 0.4 0.3333 0.99039
## Empenigr 0.004536 0.010325 0.4393 0.0 0.3333 0.99668
## Vicilath 0.002399 0.005461 0.4393 0.0 0.1667 1.00000
## Chenalbu 0.000000 0.000000    NaN 0.0 0.0000 1.00000
## Cirsarve 0.000000 0.000000    NaN 0.0 0.0000 1.00000
## Permutation: free
## Number of permutations: 0




##Q29

Question 29: Diagnostic species for a priori defined groups

Use the resulting data table to find the two top indicator species for each of the four Management types!

In this exercise we will investigate if any particular species are indicative for the four different management types. The CANOCO version of this analysis is a simplified version of these types of tests. The full versions also give significance testing and more comprehensive result presentations. See for instance IndVal by Dufresne & Legendre. 1. Go to menu Data /Add new table(s) / Indicator values. 2. Mark the Species data table, and Management as the factor defining groups. 3. Opt for Indicator Value (quantitative). 4. Export the table to Excel to be able to sort the data to find the highest values = best indicators for a group.

library(labdsv)
library(tibble)
iva <- indval(dune, dune.env.original$Management) #this time we don't need dummy variables for management
iva.df <- as.data.frame(iva$indval)
#arrange in descending order to find best indicator species, change BF to other management types as needed
arrange(rownames_to_column(iva.df), desc(HF))



##Q30

Question 30: Molecular Data

Molecular

#to follow, Stefan


##Q31

Question 31: PLS

Examine the p-values to explore which response variables actually give a significant result.

Unless you are already familiar with PLS in R we suggest that you follow the exercise in SIMCA and then try performing the same analysis by modifying the R code below. Remember you can always read the helpfile for any R function by entering ? and the function name in the console, such as ?plsr in this case.

data.pls <- readRDS("Trend_Lakes_2015_PLS.RDS")
#View(data.pls)
library(pls)
#PLS regression, like PCA, seeks to find components which maximize the variability of predictors but differs from PCA as PLS requires #the components to have maximum correlation with the response.
#Change 'Antal arter' (number of species in Swedish) to your choice of response variable(s)
pls.fit <- plsr(`Antal arter` ~ ., data = data.pls[c(9:35)], scale = TRUE, validation = "CV")
#The predictor variables are mapped to a smaller set of variables, and within that smaller space we perform a regression against the #outcome variable. PLS aims to choose new mapped variables that maximally explain the outcome variable.
#The next step is to remove unwanted variables and then build a model.  
#Cross validation is used to find the optimal number of retained dimensions.
#Then the model is rebuilt with this optimal number of dimensions. 
#Find the number of dimensions with lowest cross validation error
cv <- RMSEP(pls.fit)
best.dims <- which.min(cv$val[estimate = "adjCV", , ]) - 1
best.dims 
## 2 comps 
##       2
# Rerun the model
pls.fit2 <-  plsr(`Antal arter` ~ ., data = data.pls[c(6:35)], ncomp = best.dims)
#Finally, we extract the useful information and format the output.
coefficients <-  coef(pls.fit2)
sum.coef <-  sum(sapply(coefficients, abs))
coefficients <-  coefficients * 100 / sum.coef
coefficients <-  sort(coefficients[, 1 , 1])
barplot(tail(coefficients, 5))

#The regression coefficients are normalized so their absolute sum is 100 and the result is sorted.
#The results below show that 'Tot-N_TNb' and 'Mn' are the two most importnat positive predictors of 'Antal arter'.  You could run the code #barplot(head(coefficients, 5)) to see that at the other end of the scale for negative predictors.
summary(pls.fit2)
## Data:    X dimension: 32 29 
##  Y dimension: 32 1
## Fit method: kernelpls
## Number of components considered: 2
## TRAINING: % variance explained
##              1 comps  2 comps
## X             61.778   98.653
## Antal arter    7.295    8.371
plot(pls.fit)

validationplot(pls.fit2, val.type = "MSEP")

predplot(pls.fit2)

coefplot(pls.fit2)